{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "3.Poisson1D_Dirichlet.ipynb",
      "provenance": [],
      "collapsed_sections": [
        "Y61YA90-WIZ1",
        "ZHYqtyYJs3Jv"
      ],
      "machine_shape": "hm",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/jdtoscano94/Learning-PINNs-in-Pytorch/blob/main/3_Poisson1D_Dirichlet.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Libraries"
      ],
      "metadata": {
        "id": "Y61YA90-WIZ1"
      }
    },
    {
      "cell_type": "code",
      "source": [
        " ! pip install pyDOE"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "WZrPCr9GWEAJ",
        "outputId": "7ac89e94-8e0e-4512-e34a-68d9099da382"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Requirement already satisfied: pyDOE in /usr/local/lib/python3.7/dist-packages (0.3.8)\n",
            "Requirement already satisfied: numpy in /usr/local/lib/python3.7/dist-packages (from pyDOE) (1.21.5)\n",
            "Requirement already satisfied: scipy in /usr/local/lib/python3.7/dist-packages (from pyDOE) (1.4.1)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import torch\n",
        "import torch.autograd as autograd         # computation graph\n",
        "from torch import Tensor                  # tensor node in the computation graph\n",
        "import torch.nn as nn                     # neural networks\n",
        "import torch.optim as optim               # optimizers e.g. gradient descent, ADAM, etc.\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "import matplotlib.gridspec as gridspec\n",
        "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
        "from mpl_toolkits.mplot3d import Axes3D\n",
        "import matplotlib.ticker\n",
        "from sklearn.model_selection import train_test_split\n",
        "\n",
        "import numpy as np\n",
        "import time\n",
        "from pyDOE import lhs         #Latin Hypercube Sampling\n",
        "import scipy.io\n",
        "\n",
        "#Set default dtype to float32\n",
        "torch.set_default_dtype(torch.float)\n",
        "\n",
        "#PyTorch random number generator\n",
        "torch.manual_seed(1234)\n",
        "\n",
        "# Random number generators in other libraries\n",
        "np.random.seed(1234)\n",
        "\n",
        "# Device configuration\n",
        "device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n",
        "\n",
        "print(device)\n",
        "\n",
        "if device == 'cuda': \n",
        "    print(torch.cuda.get_device_name()) \n"
      ],
      "metadata": {
        "id": "OdYMzuSDi7Js",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "b0790762-b875-4f2c-ca4c-846a80e921ad"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "cuda\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Tunning Parameters"
      ],
      "metadata": {
        "id": "UoqCzgLSp61M"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "steps=5000\n",
        "lr=1e-3\n",
        "layers = np.array([1,50,50,20,50,50,1]) #5 hidden layers\n",
        "# To generate new data:\n",
        "min=-1\n",
        "max=1\n",
        "total_points=500\n",
        "#Nu: Number of training points (2 as we onlt have 2 boundaries), # Nf: Number of collocation points (Evaluate PDE)\n",
        "Nu=2\n",
        "Nf=250"
      ],
      "metadata": {
        "id": "L2r-CxAnp5Ub"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Problem Setup\n",
        "\n",
        "**Laplace Operator($\\Delta$)**\n",
        "\n",
        "$$\\Delta u=u_{xx}+u_{yy} = \\frac{\\partial^2 u}{\\partial x^2}+\\frac{\\partial^2 u}{\\partial y^2}$$\n",
        "\n",
        "**Harmonic Functions**\n",
        "\n",
        "$$\\Delta u =0$$\n",
        "\n",
        "They are smooth\n",
        "\n",
        "You need boundary conditions: you need the rate or the value.\n",
        "\n",
        "You have to solve all points at the same time.\n",
        "\n",
        "**Poisson equation**\n",
        "\n",
        "$$\\Delta u =f(x)$$\n",
        "\n",
        "Example: Candle in a room\n",
        "\n",
        "### Problem\n",
        "\n",
        "$$-\\Delta u=\\pi^2sin(\\pi x)$$\n",
        "\n",
        "### Boundary Conditions:\n",
        "\n",
        "$$u(-1)=0, u(1) =0$$\n",
        "\n",
        "### Exact solution:\n",
        "\n",
        "$$u(x)=sin(\\pi x)$$"
      ],
      "metadata": {
        "id": "koofMp-5Y-UA"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### PDE:\n",
        "$$-\\Delta u=\\pi^2sin(\\pi x)$$\n",
        "$$u_{xx}=-\\pi^2sin(\\pi x)$$\n",
        "$$\\frac{\\partial^2u}{\\partial x^2}=-\\pi^2sin(\\pi x)$$\n",
        "\n",
        "So the residual will be:\n",
        "\n",
        "$$0=\\frac{\\partial^2u}{\\partial x^2}+\\pi^2sin(\\pi x)$$\n",
        "\n",
        "Note: Remeber that our neural network $NN(x)\\approx u(x)$ so:\n",
        "\n",
        "$$\\frac{\\partial^2NN}{\\partial x^2}\\approx\\frac{\\partial^2u}{\\partial x^2}$$\n",
        "\n",
        "$$f=\\frac{\\partial^2NN}{\\partial x^2}+\\pi^2sin(\\pi x)\\rightarrow 0$$\n",
        "\n",
        "\n",
        "**Initial Conditions (Dirichlet BC)**\n",
        "$$u(-1)=0$$\n",
        "\n",
        "$$u(1)=0$$\n",
        "\n",
        "Set a function to describe our boundary conditions:\n",
        "\n",
        "$$f_{BC}(x)=1-|x|$$\n",
        "\n",
        "Note: $f_{BC}$ is not the answer to our PDE, but it helps us describe the boundary condition (Generate data)."
      ],
      "metadata": {
        "id": "XYwWBwmbp3S6"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Functions"
      ],
      "metadata": {
        "id": "K40IpAXPjS4I"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def f_BC(x):\n",
        "  return 1-torch.abs(x)\n",
        "def f_real(x):\n",
        "  return torch.sin(np.pi*x)\n",
        "def PDE(x):\n",
        "  return -1*(np.pi**2)*torch.sin(np.pi*x)\n"
      ],
      "metadata": {
        "id": "hZFwvWyGs2RE"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Neural Network"
      ],
      "metadata": {
        "id": "ZHYqtyYJs3Jv"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "class FCN(nn.Module):\n",
        "    ##Neural Network\n",
        "    def __init__(self,layers):\n",
        "        super().__init__() #call __init__ from parent class \n",
        "        'activation function'\n",
        "        self.activation = nn.Tanh()\n",
        "        'loss function'\n",
        "        self.loss_function = nn.MSELoss(reduction ='mean')\n",
        "        'Initialise neural network as a list using nn.Modulelist'  \n",
        "        self.linears = nn.ModuleList([nn.Linear(layers[i], layers[i+1]) for i in range(len(layers)-1)]) \n",
        "        self.iter = 0\n",
        "        'Xavier Normal Initialization'\n",
        "        # std = gain * sqrt(2/(input_dim+output_dim))\n",
        "        for i in range(len(layers)-1):\n",
        "            \n",
        "            # weights from a normal distribution with \n",
        "            # Recommended gain value for tanh = 5/3?\n",
        "            nn.init.xavier_normal_(self.linears[i].weight.data, gain=1.0)\n",
        "            \n",
        "            # set biases to zero\n",
        "            nn.init.zeros_(self.linears[i].bias.data)   \n",
        "    'foward pass'\n",
        "    def forward(self,x):\n",
        "        if torch.is_tensor(x) != True:         \n",
        "            x = torch.from_numpy(x)                \n",
        "        a = x.float()\n",
        "        for i in range(len(layers)-2):  \n",
        "            z = self.linears[i](a)              \n",
        "            a = self.activation(z)    \n",
        "        a = self.linears[-1](a)\n",
        "        return a\n",
        "    'Loss Functions'\n",
        "    #Loss BC\n",
        "    def lossBC(self,x_BC,y_BC):\n",
        "      loss_BC=self.loss_function(self.forward(x_BC),y_BC)\n",
        "      return loss_BC\n",
        "    #Loss PDE\n",
        "    def lossPDE(self,x_PDE):\n",
        "      g=x_PDE.clone()\n",
        "      g.requires_grad=True #Enable differentiation\n",
        "      f=self.forward(g)\n",
        "      f_x=autograd.grad(f,g,torch.ones([x_PDE.shape[0],1]).to(device),retain_graph=True, create_graph=True)[0] #first derivative\n",
        "      f_xx=autograd.grad(f_x,g,torch.ones([x_PDE.shape[0],1]).to(device), create_graph=True)[0]#second derivative\n",
        "      return self.loss_function(f_xx,PDE(g))\n",
        "      \n",
        "    def loss(self,x_BC,y_BC,x_PDE):\n",
        "      loss_bc=self.lossBC(x_BC,y_BC)\n",
        "      loss_pde=self.lossPDE(x_PDE)\n",
        "      return loss_bc+loss_pde"
      ],
      "metadata": {
        "id": "7EZjlbVMi8w0"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Generate data"
      ],
      "metadata": {
        "id": "pOe9yRvxjakj"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# get the analytical solution over the full domain\n",
        "x = torch.linspace(min,max,total_points).view(-1,1) #prepare to NN\n",
        "y = f_real(x)\n",
        "print(x.shape, y.shape)\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Z1V03CRRjRML",
        "outputId": "0812fa14-189d-4a77-bc88-ddc4deba672c"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "torch.Size([500, 1]) torch.Size([500, 1])\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "fig, ax1 = plt.subplots()\n",
        "ax1.plot(x.detach().numpy(),y.detach().numpy(),color='blue',label='Real Function')\n",
        "#ax1.plot(x_train.detach().numpy(),yh.detach().numpy(),color='red',label='Pred_Train')\n",
        "ax1.set_xlabel('x',color='black')\n",
        "ax1.set_ylabel('f(x)',color='black')\n",
        "ax1.tick_params(axis='y', color='black')\n",
        "ax1.legend(loc = 'upper left')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 296
        },
        "id": "OE4RUZ27p0fT",
        "outputId": "578924aa-4173-4da8-e508-3e61d081d5a2"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.legend.Legend at 0x7fb72a98c0d0>"
            ]
          },
          "metadata": {},
          "execution_count": 126
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZzV8/7A8de73ZbWS5SaCKWVudmypQhRiEpddS/6hbJEivaIpEQqSbnZriK6Ct1uWlwkmiht0qRLkyjJEu3z/v3x+Y57mmaamTPnnM9Z3s/H4zzmnO/3e873Pd85c97ns4uqYowxxhRVCd8BGGOMSUyWQIwxxoTFEogxxpiwWAIxxhgTFksgxhhjwlLKdwCxVKVKFa1Vq5bvMIwxJqEsXbr0B1Wtmnt7SiWQWrVqkZGR4TsMY4xJKCLydV7brQrLGGNMWCyBGGOMCYslEGOMMWFJqTaQvOzdu5esrCx27drlOxQToly5clSvXp3SpUv7DsUYk4+UTyBZWVkcddRR1KpVCxHxHY4BVJVt27aRlZVFWlqa73CMMfnwWoUlIs+JyBYRWZnPfhGRMSKSKSKfi8jpIfu6iMi64NYl3Bh27dpF5cqVLXnEERGhcuXKVio0Js75bgOZArQ6xP7LgDrBrRvwNICIVAIGAWcCTYFBIlIx3CAsecQf+5sYE/+8VmGp6n9EpNYhDmkDvKBuzvnFIlJBRKoBFwJzVfVHABGZi0tEr0Q3YmNMrO3eDatXw/r1sGUL/PgjiEDp0lCxItSsCbVru1sJ31+JU0y8t4EcD2wMeZwVbMtv+0FEpBuu9MIJJ5wQnSiLqWTJkjRo0IB9+/aRlpbGiy++SIUKFYr8OlOmTCEjI4OxY8cetL13794cf7y7RA0bNuSFF16ISOwADz/8MA888MAfj8855xwWLVoUsdc3qWXfPli0CN55B+bMgZUr3baCVKwIZ50FLVvCtddCnP67J5Wkz9eqOlFV01U1vWrVg0bix4XDDjuMZcuWsXLlSipVqsS4ceMifo727duzbNkyli1bFtHkAS6BhLLkYcKxYQP07+9KFBdcAKNGQYUKcN99MG0aLFsG337rSiS7d8OOHfD11/Dee/Dssy5pbNgAvXq51zj7bHjhBXesiY54TyCbgBohj6sH2/LbnvDOPvtsNm1yv8r69etp1aoVZ5xxBueddx5ffPEFALNmzeLMM8+kSZMmtGjRgu+//77I51m4cCGtW7f+43GPHj2YMmUK4KZ8GTRoEKeffjoNGjT447w7duzgr3/9Kw0aNKBhw4a8/vrr9O3bl507d9K4cWM6deoEwJFHHgm43lS9e/emfv36NGjQgGnTpv1x7gsvvJB27dpx6qmn0qlTJ2xlzNS1Zg107gwnnQSPPAKNG8Orr8K2bbBgAQwbBtdfD40aQbVqUKaMux1xhCtlnH8+3HyzSyJr1sC6dTB8OPz8M3TpAjVquMe//eb7N00+8V6FNRPoISJTcQ3mP6vqZhGZAzwc0nB+CXB/cU92113uW04kNW4MTzxRuGP379/PvHnzuOmmmwDo1q0bEyZMoE6dOnz88cfcdtttzJ8/n2bNmrF48WJEhEmTJjFixAhGjRp1yNeeNm0aH3zwAQB33nlngd1jq1Spwqeffsr48eMZOXIkkyZN4sEHH+Too49mxYoVAGzfvp1rr72WsWPHsiyPC/fGG2+wbNkyli9fzg8//MCf//xnzj//fAA+++wzVq1axXHHHce5557Lhx9+SLNmzQp3oUxS2LIFHngAnnsODjsM7r4b7rzTfeAXx0knQZ8+ruQybx48/jjcfz88+SQMGgS33AIlS0bmd0h1XhOIiLyCaxCvIiJZuJ5VpQFUdQLwDnA5kAn8Dvw12PejiDwILAleamhOg3oiyvkGv2nTJurWrUvLli3ZsWMHixYt4rrrrvvjuN1BWTwrK4v27duzefNm9uzZU6ixEu3btz+gbWThwoWHPP6aa64B4IwzzuCNN94A4N1332Xq1Kl/HFOx4qE7vn3wwQd07NiRkiVLcswxx3DBBRewZMkSypcvT9OmTalevToAjRs35r///a8lkBSRnQ3jx7vqqt9+c4mjb1+IdA2zCLRo4W4ffuiSyK23uoQ1caL7cmeKx3cvrI4F7Ffg9nz2PQc8F8l4CltSiLScNpDff/+dSy+9lHHjxtG1a1cqVKiQ5zf7nj170qtXL6666ioWLlzI4MGDi3zOUqVKkZ2d/cfj3GMuypYtC7gG/n2FacEsopzXj+Y5TPz55hvo2tVVTV1yiSsVnHpq9M977rmurWTqVFfTkJ4O/frBgAFQKt7rYeJYvLeBpJTDDz+cMWPGMGrUKA4//HDS0tJ47bXXANeesHz5cgB+/vnnP3pUPf/882Gdq2bNmqxevZrdu3fz008/MW/evAKf07JlywMa+Ldv3w5A6dKl2bt370HHn3feeUybNo39+/ezdetW/vOf/9C0adOw4jWJ7/XXoUEDWLIEJk2Cf/0rNskjhwh07AhffAGdOsHQoXDhha4h3oTHEkicadKkCQ0bNuSVV17h5ZdfZvLkyTRq1IjTTjuNN998E4DBgwdz3XXXccYZZ1ClSpWwzlOjRg2uv/566tevz/XXX0+TJk0KfE7//v3Zvn079evXp1GjRixYsABwbTUNGzb8oxE9x9VXX03Dhg1p1KgRzZs3Z8SIERx77LFhxWsS1759rj2iXTuoWxc+/xxuusl9oPtQsSI8/zy8/LKL5fTTXVuJKTpJpd4v6enpmntBqTVr1lC3bl1PEZlDsb9N4vvlF5c45s6F225zDdohtZfeZWZCmzawdi2MHg09evhLbPFMRJaqanru7VYCMcZExbffui62CxbA5Mkwblx8JQ9wPbY++giuuALuuAPuucc18pvCseYjY0zErVkDrVq5aUfeegsuvdR3RPkrXx5mzHC9wUaPduNPJk1yU6WYQ7MEgmugtsn74ksqVa0mmxUroHlzN9bivfdcG0O8K1HC9cKsUgUGDoTt293o98MO8x1ZfEv5Kqxy5cqxbds2+8CKIznrgZQrV853KKaIcpJH2bLw/vuJkTxyiLhuvePHu1JTu3Y2DUpBUr4EUr16dbKysti6davvUEyInBUJTeLISR5lyrh2jzp1fEcUnltvdWNDunWDDh3ctCpWnZW3lE8gpUuXtlXvjCmmr75ys+CWKQMLFyZu8shxyy2wa5drWP/LX1yXX5v+5GApn0CMMcWzZYtrJN+7N7FLHrn17OmSyH33ubEj48dbF9/cLIEYY8K2Ywdcfjls2uQG4yXbsJ3eveGHH2DECEhLc8nE/I8lEGNMWPbtg+uuczNY//Ofbv2NZPTII266kz593PTxHTr4jih+WAIxxoSld283n9XEiRCytEzSKVECpkxxAyO7dIHq1cEmjnZSvhuvMabonnvOjZu4807X4JzsypVzpayaNd3Kh1lZviOKD5ZAjDFF8uGH0L27W2dj5Ejf0cROpUouifz+O1xzjWtgT3WWQIwxhfbdd+4beM2abqR2qq2lUa8evPiim5L+1lsh1ccfe00gItJKRNaKSKaI9M1j/2gRWRbcvhSRn0L27Q/ZNzO2kRuTevbvhxtucDPszpjhvpGnorZt3XQnU6bAhAm+o/HL2/cHESkJjANaAlnAEhGZqaqrc45R1btDju8JhC5asVNVbVFKY2JkyBA3zuPvf4f69X1H49egQa4UcvfdcM450KiR74j88FkCaQpkqupXqroHmAq0OcTxHYFXYhKZMeYAc+fCQw+55Wi7dvUdjX8lSrhFqSpVgvbt3druqchnAjke2BjyOCvYdhARqQmkAfNDNpcTkQwRWSwibfM7iYh0C47LsPmujCm6b791S8DWq+fW9DBO1apuipMvv3QLUaWiRGlE7wBMV9X9IdtqBitk3QA8ISIn5vVEVZ2oqumqml61atVYxGpM0sjOhhtvdN+wX3sNDj/cd0Tx5aKLoH9/1x7y0ku+o4k9nwlkE1Aj5HH1YFteOpCr+kpVNwU/vwIWcmD7iDEmAp56yk1R8sQTyTdNSaQMHAjnnee6Nq9f7zua2PKZQJYAdUQkTUTK4JLEQb2pRORUoCLwUci2iiJSNrhfBTgXWJ37ucaY8K1e7abvaN0abr7ZdzTxq1QpV5VVqpRrH9q/v8CnJA1vCURV9wE9gDnAGuBVVV0lIkNF5KqQQzsAU/XAFZ/qAhkishxYAAwP7b1ljCmePXugc2c46ii3vKvNQntoNWrAmDHwwQeutJYqJJVW4ktPT9eMjAzfYRgT9/r1g4cfduM92ubbRcWEUnUj1GfPhqVL4bTTfEcUOSKyNGhzPkCiNKIbY2Jk8WIYPhz+9jdLHkUhAs8840ptXbq49VGSnSUQY8wfdu+Gm26C44+H0aN9R5N4/vQnl0SWLnUluGRnCcQY84dHHnGN5xMmQPnyvqNJTNdc46Z8GTYMVq3yHU10WQIxxgCwcqX71typk1tl0ITviSdcAr7lFjeWJllZAjHGsH+/q7o6+ujU6kUULVWruirAjz6Cp5/2HU30WAIxxjBmDHzyiRs4WKWK72iSQ+fO0LIl3H9/8i5AZQnEmBT3zTduOo7Wrd3EgCYyRFxb0r59cNttybl2iCUQY1Lc3Xe7D7exY23AYKTVrg1Dh8KsWfD6676jiTxLIMaksNmz4Y03YMAAt8qgiby77nLrhdx9N+zY4TuayLIEYkyK2rULevaEU06Be+7xHU3yKlXKTYOfleW69iYTSyDGpKhHH3Wzx44bB2XK+I4muZ17rhudPmoUrF3rO5rIsQRiTApav94NGuzQAS6+2Hc0qeHRR916KnfckTwN6pZAjEkxqu5DrEwZ943YxMYxx7gG9X//27U7JQNLIMakmLfegnfegSFD4LjjfEeTWm67DRo2dA3qybCOuiUQY1LInj2uwfzUU1N3HW+fchrUN25MjskWLYEYk0LGjYN16+Dxx6F0ad/RpKZmzdwo9VGj4OuvfUdTPF4TiIi0EpG1IpIpIn3z2N9VRLaKyLLgdnPIvi4isi64dYlt5MYknh9+cNVWrVrBZZf5jia1PfwwlCgBfQ/61Ess3hKIiJQExgGXAfWAjiJSL49Dp6lq4+A2KXhuJWAQcCbQFBgkIhVjFLoxCWnQIDeQzRrO/atRA+69F6ZOdRMuJiqfJZCmQKaqfqWqe4CpQJtCPvdSYK6q/qiq24G5QKsoxWlMwlu50s3L1L071Mvra5qJufvug2OPhV69Erdbr88EcjywMeRxVrAtt2tF5HMRmS4iNYr4XESkm4hkiEjG1q1bIxG3MQlF1X1IlS8Pgwf7jsbkOPJINzJ98WKYNs13NOGJ90b0WUAtVW2IK2U8X9QXUNWJqpququlVq1aNeIDGxLt33oG5c10Vlk3VHl+6dIHGjaFPH9i503c0ReczgWwCaoQ8rh5s+4OqblPV3cHDScAZhX2uMQb27nXddk8+2Y1BMPGlZEnXJvXNN/Dkk76jKTqfCWQJUEdE0kSkDNABmBl6gIhUC3l4FbAmuD8HuEREKgaN55cE24wxISZPdnMvPfaYzXcVr5o3h6uucj2zvv/edzRF4y2BqOo+oAfug38N8KqqrhKRoSJyVXDYHSKySkSWA3cAXYPn/gg8iEtCS4ChwTZjTGDHDtfm0awZXHml72jMoYwY4aqwBg3yHUnRiCZq838Y0tPTNSMjw3cYxsTEQw+5dT4+/BDOOcd3NKYgPXq4nnKrV7sqx3giIktVNT339nhvRDfGhGHrVvettm1bSx6JYsAAKFfOLS+cKCyBGJOEhg1zk/Ulw3xLqeKYY1x369deg0SpKLEEYkyS2bABxo+Hv/0N6tb1HY0pinvvdV2t+/RJjMGFlkCMSTIDBrjuoTZoMPGUL++qsObPd2N34p0lEGOSyGefwcsvw113wfF5zs1g4l337lCrlptoMTvbdzSHZgnEmCRy//1QsaKrAjGJqWxZePBB92Xg1Vd9R3NolkCMSRLz58OcOdCvH1So4DsaUxw33OBWLuzf3y0CFq8sgRiTBFThgQegenW4/Xbf0ZjiKlECHnkE1q+HSZN8R5M/SyDGJIF33oGPP4aBA91YApP4LrsMzjvPdcmO14kWLYEYk+Cys13Pq9q1oWtX39GYSBFxbSHffutGqMcjSyDGJLgZM1yD66BBts55srngAmjRwlVn7djhO5qDWQIxJoHt3++qrU49FTp18h2NiYYHH3RT04wd6zuSg1kCMSaBTZvmJt8bMsQNHjTJ56yz4Ior3NxmP//sO5oDWQIxJkHt2+eqrRo2hHbtfEdjomnoUNi+HZ54wnckB7IEYkyCeuEFyMx0Hy4l7D85qZ1+OlxzDTz+OPwYRysfeX3biUgrEVkrIpki0jeP/b1EZLWIfC4i80SkZsi+/SKyLLjNzP1cY5LZnj0ucaSnu9XsTPIbMgR+/RVGjvQdyf94SyAiUhIYB1wG1AM6iki9XId9BqSrakNgOjAiZN9OVW0c3OxfyKSUyZPh66/dolEivqMxsVC/PnTo4NZO37LFdzSOzxJIUyBTVb9S1T3AVKBN6AGqukBVfw8eLgaqxzhGY+LOzp0ucTRrBpdc4jsaE0uDBsGuXfDoo74jcXwmkOOBjSGPs4Jt+bkJmB3yuJyIZIjIYhFpm9+TRKRbcFzG1q1bixexMXFgwgQ3uMxKH6nnlFPgxhvdei/ffus7mgRpRBeRzkA68FjI5prBGr03AE+IyIl5PVdVJ6pquqqmV61aNQbRGhM9v//uvn02b+4GmZnUM3Cg64EXD6tN+kwgm4AaIY+rB9sOICItgH7AVaq6O2e7qm4Kfn4FLASaRDNYY+LBs8/C99+7qgyTmtLS3JQ1zz4Lmw76xIwtnwlkCVBHRNJEpAzQATigN5WINAGewSWPLSHbK4pI2eB+FeBcYHXMIjfGg5y67wsvhPPP9x2N8emBB9wcaL7bQrwlEFXdB/QA5gBrgFdVdZWIDBWRnF5VjwFHAq/l6q5bF8gQkeXAAmC4qloCMUlt0iTYvNlVYZjUlpbm2kImTvTbFiKaCCu3R0h6erpmZGT4DsOYItu9G0480c24+9571nhu4Kuv4OSToUeP6I9QF5GlQZvzARKiEd2YVPfcc66+e+BASx7GqV3blUKeecaVTH2wBGJMnNu9203nfc45cPHFvqMx8aRfP9i710206IMlEGPi3JQpsHGj63llpQ8T6sQToXNnNzbou+9if35LIMbEsT17XH//M8+Eli19R2PiUb9+7n3y2GMFHxtplkCMiWMvvADffGOlD5O/OnXcYmJPP+3GCMWSJRBj4tTeva70kZ4OrVr5jsbEs/79XVtZrGfqtQRiTJx66SXYsMFKH6ZgJ58MN9zg5siK5Uy9lkCMiUP79sGwYW4hoSuu8B2NSQT9+7vZCkaNit05LYEYE4f+8Q9Yv97GfZjCO+UUt17IuHEQq4nHLYEYE2f273dTtTdqZKsNmqLp39/N2ByrUoglEGPizNSpsG6dlT5M0dWtC+3bw9ix8MMP0T+fJRBj4sj+/fDgg9CgAbTNd5k0Y/I3YIArhTz+ePTPZQnEmDjy2muwdq37EChh/50mDPXqwXXXuVLIjz9G91z2FjUmTmRnu9JHvXpw7bW+ozGJrH9/+PVXePLJ6J6nVGEOEpE/4RZtOg7YCawEMlQ1O4qxGZNSXn8dVq+GV16x0ocpngYN4OqrXQLp1QuOPjo65znk21RELhKROcDbwGVANaAe0B9YISJDRKR8dEIzJnVkZ8PQoXDqqa76wZjiGjAAfv4Znnoqeuco6HvO5cAtqvpnVe2mqv1V9V5VvQpoBHwGhD3Fm4i0EpG1IpIpIn3z2F9WRKYF+z8WkVoh++4Ptq8VkUvDjcGYePDPf8LKle6fvmRJ39GYZNCkCbRuDaNHu+qsaDhkAlHV3qr6TT779qnqP1X19XBOLCIlgXG4kk09oKOI1Mt12E3AdlU9CRgNPBo8tx5uDfXTgFbA+OD1jEk4OaWPk092XTCNiZQBA1xD+vjx0Xn9QtW0isiLInJ0yONaIjKvmOduCmSq6lequgeYCrTJdUwb4Png/nTgYhGRYPtUVd2tqhuAzOD1ouLRR6HvQeUjYyJj1ixYvtw1fFrpw0RS06ZuIs6RI+G33yL/+oVtqvsA+FhELheRW4B/A8Vdhfd4YGPI46xgW57HqOo+4GegciGfC4CIdBORDBHJ2Brm+P7//tf1qd64scBDjSkSVRgyBE46CTp29B2NSUYDBkCFCm5izkgrVAJR1WeAm4E3gaHA+ao6K/LhRJ6qTlTVdFVNr1q1alivcf/97ufw4REMzBjgrbfgs8/cokClCtUn0piiOecc+OILqF8/8q9d2CqsvwDPATcCU4B3RKRRMc+9CagR8rh6sC3PY0SkFHA0sK2Qz42YE06Arl1h0iTIyorWWUyqUXVtH2lpbkEgY6IlWlWjha3CuhZopqqvqOr9QHdcIimOJUAdEUkTkTK4RvGZuY6ZCXQJ7rcD5quqBts7BL200oA6wCfFjOeQHnjANXb6WrzeJJ/ZsyEjw5U+Spf2HY0xRVfYKqy2qrol5PEnwJnFOXHQptEDmAOsAV5V1VUiMlREcuYgnQxUFpFMoBfQN3juKuBVYDXwL+B2Vd1fnHgKUqsWdOkCEyfCt99G80wmFeSUPmrWhBtv9B2NMeER94U+n50i/YHxqprnjCoi0hw4XFXfilJ8EZWenq4ZGRlhP/+rr1xXyx494InidiEwKW3OHNc75plnoFs339EYc2gislRV03NvL6jZbgUwS0R2AZ8CW4FyuCqjxsC7wMMRjjVu1a4Nf/mL+6fv2xeOPdZ3RCYR5fS8ymlbMyZRFVSF1U5Vz8VVM60CSgK/AC8BTVX1blWN0dpX8aFfP9izBx57zHckJlHNmwcffeR695Up4zsaY8JXUBXWaqAFMBu4KPf+/Kq24lVxq7By3HgjTJ/u+lUfc0wEAjMpQxXOP9+NLcrMhLJlfUdkTMHyq8IqqAQyAZgHnApkhNyWBj9TUv/+sHt3bBevN8lhwQL44ANXBWrJwyS6Q5ZA/jhI5GlVvTUG8URVpEogAJ07w4wZ7ptkmOMTTQq68EK3XO369VCunO9ojCmccEsgACRD8oi0fv1g504rhZjCe+89d+vTx5KHSQ62bE2YYr14vUl8Q4a4nnu33OI7EmMiwxJIMfTv7xavHz3adyQm3r3/vmv/6NMHDjvMdzTGRIYlkGI47TRo186t+BXtxetNYhs61PXYs0GDJplYAimmAQPcal82Mt3kZ9EiePdd6N0bDj/cdzTGRI4lkGJq0ACuvdYtXr99u+9oTDwaMsT11Ove3XckxkSWJZAIGDAAfvkFxozxHYmJN4sXw7//7UofRxzhOxpjIssSSAQ0agRt27pqrJ9/9h2NiSdDh0KVKnCrdYQ3ScgSSIQMHAg//WSlEPM/n3zi1vy45x448kjf0RgTeZZAIqRJE7jyStel95dffEdj4sGDD0KlSnD77b4jMSY6LIFE0MCBriF97FjfkRjfMjLceue9esFRR/mOxpjo8JJARKSSiMwVkXXBz4p5HNNYRD4SkVUi8rmItA/ZN0VENojIsuDWOLa/Qd7S0+Hyy930Jr/+6jsa49PAga700bOn70iMiR5fJZC+wDxVrYOb7bdvHsf8DtyoqqcBrYAnRKRCyP7eqto4uC2LfsiFM2iQG1T41FO+IzG+fPSRa/u47z4oX953NMZEj68E0gZ4Prj/PNA29wGq+qWqrgvufwtsAeJ+3tumTeGKK2DkSOuRlaoGDnTjPqztwyQ7XwnkGFXdHNz/Djjkskwi0hQoA6wP2TwsqNoaLSL5rqwgIt1EJENEMrZujc3iiUOHurYQmyMr9fznP27Ued++1vPKJL9CrQcS1guLvAvktWp4P+B5Va0Qcux2VT2oHSTYVw1YCHRR1cUh277DJZWJwHpVHVpQTJFcD6Qg7dq5AWQbNkDlyjE5pfFM9cD1PmzSRJMsirUeSDhUtYWq1s/j9ibwfZAEcpLBlnyCLg+8DfTLSR7Ba29WZzfwd6BptH6PcA0ZAjt22NrpqWTePFcCeeABSx4mNfiqwpoJdAnudwHezH2AiJQBZgAvqOr0XPtyko/g2k9WRjXaMJx2GtxwgxtY+N13vqMx0abqprSpXt3W+zCpw1cCGQ60FJF1QIvgMSKSLiKTgmOuB84HuubRXfdlEVkBrACqAA/FNvzCGTQI9uyB4cN9R2KibfZsN+9V//621rlJHVFrA4lHsWwDyXHzzfDSS5CZ6b6dmuSjCn/+M2zbBmvXQpkyviMyJrJi3gZinAEDIDsbhg3zHYmJlpkzYelS133XkodJJZZAoqxmTVcnPmmS65Flkkt2tkscJ50Ef/mL72iMiS1LIDHQrx+UKuXGh5jkMn06fP65a+8qVcp3NMbEliWQGDjuOLjtNnjhBVdHbpLD3r2u0bx+fejY0Xc0xsSeJZAY6dPHjQ0YPNh3JCZSnnvODRp8+GEoWdJ3NMbEniWQGPnTn+DOO2HaNFi+3Hc0prh+/90NFj33XGjd2nc0xvhhCSSGeveGChXcPEkmsY0ZA5s3uzE+Ir6jMcYPSyAxVKGCa1D/179g/nzf0Zhwbd8Ojz7qZl1u1sx3NMb4Ywkkxm6/HU44wbWJZGf7jsaEY/hwN1X/ww/7jsQYvyyBxFi5cm6t7IwM1wXUJJZNm1z1VadO0LCh72iM8csSiAedOkGDBm7W1j17fEdjimLoUNi/38b0GAOWQLwoWdJVg6xfD88+6zsaU1hr18LkydC9O6Sl+Y7GGP8sgXhy2WVu8aEhQ+DXX31HYwqjf39XBdm/v+9IjIkPlkA8EXE9ebZuhVGjfEdjCrJokWuz6t3bjekxxlgC8appU7f07ciRtuhUPMvOhrvvdlPS3Huv72iMiR9eEoiIVBKRuSKyLviZ33ro+0MWk5oZsj1NRD4WkUwRmRasXpiQHn4Ydu1yVVkmPk2bBp984v5WRxzhOxpj4oevEkhfYJ6q1gHmBY/zslNVGwe3q0K2PwqMVtWTgO3ATdENN3rq1HETLU6cCCtW+I7G5LZzp5s5oEkTm/I9FiMAABGeSURBVK7dmNx8JZA2wPPB/edx65oXSrAOenMgZxRFkZ4fjwYNgqOPdtUkKbRAZEJ44gn45hvXTlXCKnyNOYCvf4ljVHVzcP874Jh8jisnIhkislhEcpJEZeAnVd0XPM4Cjs/vRCLSLXiNjK1bt0Yk+EirXNlVYc2bB7Nm+Y7G5Pj+e1dt1aYNXHSR72iMiT9RSyAi8q6IrMzj1ib0OHWLsuf3vbtmsA7vDcATInJiUeNQ1Ymqmq6q6VWrVi36LxIj3btD3bpwzz02uDBeDBzo2qdGjPAdiTHxKWoJRFVbqGr9PG5vAt+LSDWA4OeWfF5jU/DzK2Ah0ATYBlQQkZz136oDm6L1e8RK6dLw+OOQmQljx/qOxqxc6ZYhvv12OPlk39EYE598VWHNBLoE97sAb+Y+QEQqikjZ4H4V4FxgdVBiWQC0O9TzE1GrVm6A4dChbnyI8UPVtUeVL+9KIcaYvPlKIMOBliKyDmgRPEZE0kVkUnBMXSBDRJbjEsZwVV0d7OsD9BKRTFybyOSYRh9Fjz8OO3bYB5dPb7wB777rJr2sVMl3NMbEL9EU6vaTnp6uGRkZvsMo0J13umqszz6zGV9j7bffXFtUxYqwdCmUKlXwc4xJdiKyNGiPPoB1TIxDgwa5D7AePaxbb6w98ghs3AjjxlnyMKYglkDiUKVKbp6s99+HF1/0HU3qWLcOHnsMOne2lQaNKQxLIHHqr3+Fs892cy9t3+47muSn6qoOy5a1brvGFJYlkDhVogQ8/TRs2+bWUTfRNWsWzJ4NgwdDtWq+ozEmMVgCiWONGkHPnjBhAixZ4jua5PXbb670Ua+eu97GmMKxBBLnhg6FY4+FW291S6mayBsyBP77Xxg/3g3oNMYUjiWQOFe+vBsbsnSpK4mYyPrsM3d9b74ZLrjAdzTGJBZLIAmgfXto0QLuvx+ysnxHkzz274dbboEqVazh3JhwWAJJACLwzDPuA697dxsbEilPPeVKdk8+6cbdGGOKxhJIgqhdG4YNg7ffhn/8w3c0ie/rr6F/f7j8crj+et/RGJOYLIEkkJ493diQO++ELXnOX2wKQxX+7//c/fHjXQnPGFN0lkASSMmSborxX3+FO+7wHU3ievZZmDPHtXvUrOk7GmMSlyWQBFOvHgwYANOmwYwZvqNJPBs2QK9ecPHFrj3JGBM+SyAJqE8fOP106NYNvvvOdzSJIzvbTRFTogQ895ytcW5Mcdm/UAIqXRpeesmtG3LTTdYrq7Ceegree8/1ujrhBN/RGJP4LIEkqLp13cyx77zjuviaQ/viCzeOpnVr6NrVdzTGJAcvCUREKonIXBFZF/w8qBe+iFwkIstCbrtEpG2wb4qIbAjZ1zj2v4V/t98Ol14K99wDX37pO5r4tWuXG4x5xBEwcaL1ujImUnyVQPoC81S1DjAveHwAVV2gqo1VtTHQHPgd+HfIIb1z9qvqsphEHWdEXF1+uXJuDYs9e3xHFJ9694bPP4cpU2ymXWMiyVcCaQM8H9x/HmhbwPHtgNmq+ntUo0pAxx3nuqUuWQJ9D0rD5s033fLAd98NV1zhOxpjkouvBHKMqm4O7n8HHFPA8R2AV3JtGyYin4vIaBEpm98TRaSbiGSISMbWrVuLEXL8uuYaN8hw9Gh44w3f0cSPjRtdr6szznBL1RpjIks0Sl14RORd4Ng8dvUDnlfVCiHHblfVPGcjEpFqwOfAcaq6N2Tbd0AZYCKwXlWHFhRTenq6ZmRkFPl3SQR79sB557nG4k8/hRNP9B2RX3v3QvPmsGyZm3H3pJN8R2RM4hKRpaqannt7qWidUFVbHCKY70WkmqpuDpLBoSbmuB6YkZM8gtfOKb3sFpG/A/dGJOgEVqaMG1x4+ulw3XWwaJFrG0lV994LH3zg5g2z5GFMdPiqwpoJdAnudwHePMSxHclVfRUkHUREcO0nK6MQY8KpVQteeMF94+7ZM3XHh7z4IowZ49o9Onb0HY0xyctXAhkOtBSRdUCL4DEiki4ik3IOEpFaQA3gvVzPf1lEVgArgCrAQzGIOSG0bg0PPODmzBo71nc0sffpp26E/oUX2hofxkRb1NpA4lEyt4GEys52DeuzZsHs2XDJJb4jio3vvoMzz3Qlr6VLoWpV3xEZkxzyawOxkehJqEQJV41z2mlurYu1a31HFH2//QZXXgk//OAmmbTkYUz0WQJJUkcdBTNnusb1K66A77/3HVH07N8PN9zgqq+mTnXddo0x0WcJJInVquWSyObNcNll8MsvviOKPFU3PfvMmW6SxCuv9B2RManDEkiSO+ssmD7dTeVx9dWwe7fviCJr6ND/9bjq0cN3NMakFksgKeCyy9ycWfPnQ6dObpBdMhg5EgYPdrPrjhzpOxpjUo8lkBRx441uqpPXX3djIxI9iTz9tJsk8frrXZdlWxzKmNizf7sUctdd8PjjLom0b5+4s/eOGQO33ebaO156ya0Vb4yJPUsgKebuu11j84wZcO218HsCzW+sCkOGwJ13uvac115zqzMaY/ywBJKC7rjDVQG9/TZcdBFsOdRMZHEiO9slv5w2j1dfhbL5zsFsjIkFSyApqnt3N/X7ihVw9tnxvaLhr7+6kfVPPumSyOTJUCpq04AaYwrLEkgKa9sWFixwH9BnneXWV483GzbAuefCW2/BU0/BqFHWYG5MvLB/xRR35pmweDHUrOlGrA8Y4EZ2x4Pp06FJE7cw1OzZbpyHrWduTPywBGKoXdutH/K3v8FDD0GzZn6rtH79Ff7v/9y6Jqec4qYoadnSXzzGmLxZAjEAHHaYa1t4+WU3+WKjRm5wXqy7+s6aBfXquXXee/d2i0KlpcU2BmNM4VgCMQe44QZYtcp94+/dGxo0cL21oj3r/4oVbi2Tq66CChVciWjECOuma0w8swRiDlKtmpuc8O23XZtD69aup9bMma47bSR99hl07uxKPB9+CI8+6tbyOOusyJ7HGBN5XhKIiFwnIqtEJFtEDlqkJOS4ViKyVkQyRaRvyPY0Efk42D5NRMrEJvLUcvnlbhLGCRPcWJE2bVybxLBh8M034b/uTz/BlClu1cDTT4d//hPuuQfWr4f77nNT0Btj4p+XFQlFpC6QDTwD3KuqBy0TKCIlgS+BlkAWsAToqKqrReRV4A1VnSoiE4Dlqvp0QedNlRUJo2HfPjd479lnYeFCt61hQ7fa4VlnuaquE088eFqR7GzYtMm1qyxaBO+/D++95+biSktzU5LcfLOrtjLGxKf8ViT0uqStiCwk/wRyNjBYVS8NHt8f7BoObAWOVdV9uY87FEsgkbFhg5tGZM4c18id09BeooRLBBUrusSxaxf8+OP/ppAXcUmnRQs3F1d6unXLNSYR5JdA4nk87/HAxpDHWcCZQGXgJ1XdF7L9+PxeRES6Ad0ATjjhhOhEmmLS0lxV0333wc6dsHq1awRfvx62b3e3kiWhXDmXUE480d3S062kYUwyiVoCEZF3gWPz2NVPVd+M1nlzU9WJwERwJZBYnTdVHHaYW0LWlpE1JvVELYGoaotivsQmoEbI4+rBtm1ABREpFZRCcrYbY4yJoXjuxrsEqBP0uCoDdABmqmu0WQC0C47rAsSsRGOMMcbx1Y33ahHJAs4G3haROcH240TkHYCgdNEDmAOsAV5V1VXBS/QBeolIJq5NZHKsfwdjjEl1XnthxZr1wjLGmKLLrxdWPFdhGWOMiWOWQIwxxoTFEogxxpiwWAIxxhgTlpRqRBeRrcDXYT69CvBDBMOJFIuraCyuorG4iiZZ46qpqlVzb0ypBFIcIpKRVy8E3yyuorG4isbiKppUi8uqsIwxxoTFEogxxpiwWAIpvIm+A8iHxVU0FlfRWFxFk1JxWRuIMcaYsFgJxBhjTFgsgRhjjAmLJZAQInKdiKwSkWwRybfLm4i0EpG1IpIpIn1DtqeJyMfB9mnBNPSRiKuSiMwVkXXBz4p5HHORiCwLue0SkbbBvikisiFkX+NYxRUctz/k3DNDtvu8Xo1F5KPg7/25iLQP2RfR65Xf+yVkf9ng988MrketkH33B9vXikiByzZHOK5eIrI6uD7zRKRmyL48/6YxiquriGwNOf/NIfu6BH/3dSLSJcZxjQ6J6UsR+SlkX1Sul4g8JyJbRGRlPvtFRMYEMX8uIqeH7Cv+tVJVuwU3oC5wCrAQSM/nmJLAeqA2UAZYDtQL9r0KdAjuTwBujVBcI4C+wf2+wKMFHF8J+BE4PHg8BWgXhetVqLiAHfls93a9gJOBOsH944DNQIVIX69DvV9CjrkNmBDc7wBMC+7XC44vC6QFr1MyhnFdFPIeujUnrkP9TWMUV1dgbB7PrQR8FfysGNyvGKu4ch3fE3guBtfrfOB0YGU++y8HZgMCnAV8HMlrZSWQEKq6RlXXFnBYUyBTVb9S1T3AVKCNiAjQHJgeHPc80DZCobUJXq+wr9sOmK2qv0fo/Pkpalx/8H29VPVLVV0X3P8W2AIcNNI2AvJ8vxwi3unAxcH1aQNMVdXdqroByAxeLyZxqeqCkPfQYtzqn9FWmOuVn0uBuar6o6puB+YCrTzF1RF4JULnzpeq/gf3ZTE/bYAX1FmMW821GhG6VpZAiu54YGPI46xgW2XgJ3ULYYVuj4RjVHVzcP874JgCju/AwW/eYUERdrSIlI1xXOVEJENEFudUqxFH10tEmuK+Va4P2Ryp65Xf+yXPY4Lr8TPu+hTmudGMK9RNuG+yOfL6m8YyrmuDv890EclZ+jourldQ1ZcGzA/ZHK3rVZD84o7ItYramujxSkTeBY7NY1c/VfW2NO6h4gp9oKoqIvn2vQ6+XTTAreSY437cB2kZXH/wPsDQGMZVU1U3iUhtYL6IrMB9SIYtwtfrRaCLqmYHm8O+XslIRDoD6cAFIZsP+puq6vq8XyHiZgGvqOpuEfk/XOmteYzOXRgdgOmquj9km8/rFTUpl0BUtUUxX2ITUCPkcfVg2zZc8bBU8C0yZ3ux4xKR70WkmqpuDj7wthzipa4HZqjq3pDXzvk2vltE/g7cG8u4VHVT8PMrEVkINAFex/P1EpHywNu4Lw+LQ1477OuVh/zeL3kdkyUipYCjce+nwjw3mnEhIi1wSfkCVd2dsz2fv2kkPhALjEtVt4U8nIRr88p57oW5nrswAjEVKq4QHYDbQzdE8XoVJL+4I3KtrAqr6JYAdcT1ICqDe7PMVNcytQDX/gDQBYhUiWZm8HqFed2D6l6DD9Gcdoe2QJ49NqIRl4hUzKkCEpEqwLnAat/XK/jbzcDVD0/PtS+S1yvP98sh4m0HzA+uz0ygg7heWmlAHeCTYsRSpLhEpAnwDHCVqm4J2Z7n3zSGcVULeXgVsCa4Pwe4JIivInAJB5bEoxpXENupuEbpj0K2RfN6FWQmcGPQG+ss4OfgC1JkrlU0egYk6g24GlcXuBv4HpgTbD8OeCfkuMuBL3HfIPqFbK+N+wfPBF4DykYorsrAPGAd8C5QKdieDkwKOa4W7ptFiVzPnw+swH0QvgQcGau4gHOCcy8Pft4UD9cL6AzsBZaF3BpH43rl9X7BVYldFdwvF/z+mcH1qB3y3H7B89YCl0X4/V5QXO8G/wc512dmQX/TGMX1CLAqOP8C4NSQ5/4tuI6ZwF9jGVfweDAwPNfzona9cF8WNwfv5SxcW1V3oHuwX4BxQcwrCOldGolrZVOZGGOMCYtVYRljjAmLJRBjjDFhsQRijDEmLJZAjDHGhMUSiDHGmLBYAjHGGBMWSyDGGGPCYgnEGI9E5M/BpIDlROQIceuT1PcdlzGFYQMJjfFMRB7CjUY/DMhS1Uc8h2RMoVgCMcazYG6lJcAu4Bw9cBZXY+KWVWEZ419l4EjgKFxJxJiEYCUQYzwTt0b2VNwiRNVUtYfnkIwplJRbD8SYeCIiNwJ7VfUfIlISWCQizVV1fkHPNcY3K4EYY4wJi7WBGGOMCYslEGOMMWGxBGKMMSYslkCMMcaExRKIMcaYsFgCMcYYExZLIMYYY8Ly/2QaPIQGwwWsAAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "#def get_training_data(x):\n",
        "#Nu: Number of training point, # Nf: Number of colloction points\n",
        "# Set Boundary conditions x=min & x= max\n",
        "BC_1=x[0,:]\n",
        "BC_2=x[-1,:]\n",
        "# Total Tpaining points BC1+BC2\n",
        "all_train=torch.vstack([BC_1,BC_2])\n",
        "#Select Nu points\n",
        "idx = np.random.choice(all_train.shape[0], Nu, replace=False) \n",
        "x_BC=all_train[idx]\n",
        "#Select Nf points\n",
        "# Latin Hypercube sampling for collocation points \n",
        "x_PDE = BC_1 + (BC_2-BC_1)*lhs(1,Nf)\n",
        "x_PDE = torch.vstack((x_PDE,x_BC)) "
      ],
      "metadata": {
        "id": "_JMjSMOKdlWt"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "y_BC=f_BC(x_BC).to(device)"
      ],
      "metadata": {
        "id": "SxRNKJ4Ejif8"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Train Neural Network"
      ],
      "metadata": {
        "id": "I82eH55ElVSq"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#Store tensors to GPU\n",
        "torch.manual_seed(123)\n",
        "x_PDE=x_PDE.float().to(device)\n",
        "x_BC=x_BC.to(device)\n",
        "#Create Model\n",
        "model = FCN(layers)\n",
        "print(model)\n",
        "model.to(device)\n",
        "params = list(model.parameters())\n",
        "optimizer = torch.optim.Adam(model.parameters(),lr=lr,amsgrad=False)\n",
        "start_time = time.time()"
      ],
      "metadata": {
        "id": "8KZvgKcalOVV",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "fd82117b-8c0c-443b-c567-7dedd5dc2170"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "FCN(\n",
            "  (activation): Tanh()\n",
            "  (loss_function): MSELoss()\n",
            "  (linears): ModuleList(\n",
            "    (0): Linear(in_features=1, out_features=50, bias=True)\n",
            "    (1): Linear(in_features=50, out_features=50, bias=True)\n",
            "    (2): Linear(in_features=50, out_features=20, bias=True)\n",
            "    (3): Linear(in_features=20, out_features=50, bias=True)\n",
            "    (4): Linear(in_features=50, out_features=50, bias=True)\n",
            "    (5): Linear(in_features=50, out_features=1, bias=True)\n",
            "  )\n",
            ")\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "model.lossBC(x_BC,y_BC)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "FGZzo1X40SJ_",
        "outputId": "a80c6cc8-1dab-48ce-fa7d-1ff214a90781"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "tensor(0.0160, device='cuda:0', grad_fn=<MseLossBackward0>)"
            ]
          },
          "metadata": {},
          "execution_count": 119
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "for i in range(steps):\n",
        "    yh = model(x_PDE)\n",
        "    loss = model.loss(x_BC,y_BC,x_PDE)# use mean squared error\n",
        "    optimizer.zero_grad()\n",
        "    loss.backward()\n",
        "    optimizer.step()\n",
        "    if i%(steps/10)==0:\n",
        "      print(loss)"
      ],
      "metadata": {
        "id": "KM3GHnh9ljCh",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "0b35909d-eda4-4c37-8a44-83f26782ce58"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "tensor(47.8950, device='cuda:0', grad_fn=<AddBackward0>)\n",
            "tensor(0.0152, device='cuda:0', grad_fn=<AddBackward0>)\n",
            "tensor(0.0030, device='cuda:0', grad_fn=<AddBackward0>)\n",
            "tensor(0.0016, device='cuda:0', grad_fn=<AddBackward0>)\n",
            "tensor(0.0008, device='cuda:0', grad_fn=<AddBackward0>)\n",
            "tensor(0.0004, device='cuda:0', grad_fn=<AddBackward0>)\n",
            "tensor(0.0002, device='cuda:0', grad_fn=<AddBackward0>)\n",
            "tensor(0.0003, device='cuda:0', grad_fn=<AddBackward0>)\n",
            "tensor(0.0032, device='cuda:0', grad_fn=<AddBackward0>)\n",
            "tensor(0.0005, device='cuda:0', grad_fn=<AddBackward0>)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Plots"
      ],
      "metadata": {
        "id": "kpRh89zM8Te8"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Function\n",
        "yh=model(x.to(device))\n",
        "y=f_real(x)\n",
        "#Error\n",
        "print(model.lossBC(x.to(device),f_real(x).to(device)))"
      ],
      "metadata": {
        "id": "6kJ3gQ6WYhp9",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "9f07415b-9aaf-4c63-d591-ce94b6e1daa4"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "tensor(8.2365e-06, device='cuda:0', grad_fn=<MseLossBackward0>)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Derivative\n",
        "g=x.to(device)\n",
        "g=g.clone()\n",
        "g.requires_grad=True #Enable differentiation\n",
        "f=model(g)\n",
        "f_x=autograd.grad(f,g,torch.ones([g.shape[0],1]).to(device),retain_graph=True, create_graph=True)[0]"
      ],
      "metadata": {
        "id": "jEUUqEjr8Sim"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Detach from GPU\n",
        "y_plot=y.detach().numpy()\n",
        "yh_plot=yh.detach().cpu().numpy()\n",
        "f_x_plot=f_x.detach().cpu().numpy()"
      ],
      "metadata": {
        "id": "udWgFd6mqY7B"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Plot\n",
        "fig, ax1 = plt.subplots()\n",
        "ax1.plot(x,y_plot,color='blue',label='Real u(x)')\n",
        "ax1.plot(x,yh_plot,color='red',label='Predicted u(x)')\n",
        "ax1.plot(x,f_x_plot,color='green',label='Computed u\\'(x)')\n",
        "ax1.set_xlabel('x',color='black')\n",
        "ax1.set_ylabel('f(x)',color='black')\n",
        "ax1.tick_params(axis='y', color='black')\n",
        "ax1.legend(loc = 'upper left')"
      ],
      "metadata": {
        "id": "M9YZDGcpodnQ",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 296
        },
        "outputId": "61f96fb1-4c77-42dd-f094-dbeb55a9c530"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.legend.Legend at 0x7fb72a152810>"
            ]
          },
          "metadata": {},
          "execution_count": 130
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd1QU1/vH8fdQFLEAdqwIigqIKIgVsHdjLLF+LbFHo4nRJGoilthbTNBYYk80MfbYYsECFlRQxIINCxaUZkMUWfb+/jDyM4mKBbi7cF/n7Im7zM580JN55k55riaEQFEURcl+TGQHUBRFUeRQBUBRFCWbUgVAURQlm1IFQFEUJZtSBUBRFCWbMpMd4G0ULFhQ2NnZyY6hKIpiVEJCQmKFEIX+/blRFQA7OzuCg4Nlx1AURTEqmqZde9nn6hSQoihKNqUKgKIoSjalCoCiKEo2ZVTXAF4mOTmZGzdu8OTJE9lRlHdgYWFBiRIlMDc3lx1FUbIdoy8AN27cIG/evNjZ2aFpmuw4ylsQQhAXF8eNGzcoU6aM7DiKku0Y/SmgJ0+eUKBAAbXzN0KaplGgQAE1elMUSYy+AABq52/E1L+doshj9KeAFCWzJemSOBNzhst3LxP9KJr4x/FoaJiZmFHAsgAl85XE3sYeh/wOmGhZ4hhLyaJUAVCUNOj0Og5dP8TWC1vZeXknp6NPo9Pr0vyejYUN1UtUp4lDE9pUaENp69KZkFZR3pwqAOnA1NSUSpUqodPpKFOmDL/88gvW1tZvvZ5ly5YRHBzMnDlz3vg7jx8/pmnTpuzZswdTU9OXLnPq1ClmzpzJsmXL3jpTdnb13lUWHV/E0tCl3Hp4CzMTM7xKefFlrS+pUrQKjgUcKZKnCAVyFUAg0Ol1xDyKIfJ+JOfjznPkxhEOXj/I0B1DGbpjKNWLV+cTj0/o6NIRCzML2b+eosgrAJqmWQABQM6/c6wVQoyRled95MqVi9DQUAB69OjB3Llz+eabbzJl20uWLKFt27av3PkDVKpUiRs3bhAZGUmpUqUyJZcxOxd7jgkBE/jt9G8ANCvbjNlNZtOkbBPy5cz3yu/lMM1BaevSlLYujVdpL/pU7QPApfhLrA9fz7LQZfTc1JNhO4cxrOYwhlQfQu4cuTPld1KUl5E5AkgC6gshEjRNMwcOaJq2XQgR9K4r/Pxz+Hs/nG7c3GD27DdfvmbNmoSFhQEQERHBoEGDiImJwdLSkp9//pkKFSqwefNmJkyYwNOnTylQoAArV66kSJEir1zn2LFjyZMnD8OHDwfAxcWFLVu2YGdnx8qVK1m1ahUAGzZsYM6cOezevZvbt2/j4+NDQEAARYsWpVWrVvz+++989dVX7/6XkcXFPIphpP9IlpxYQi7zXHxR4wuGVB9CSauS77XesvnL8lXtr/iy1pfsvbqXWYdnMWrPKGYfmY2vty8DPAZgavLqAq4oGUXaFSrxTMLfb83/fhn1BMUpKSn4+/vzwQcfANCvXz/8/PwICQlhxowZDBw4EIA6deoQFBTEiRMn6NSpE9OmTXun7T19+pTLly/zvENqmzZtsLW1Ze7cufTt25dx48ZRtGhRADw8PAgMDHz/XzIL0gs9847Nw3GOI8tPLufzGp9z5bMrTG88/b13/i/SNI36ZeqzpcsWDvU6hFMhJz7d/inVF1Un5FZIum1HUd6U1GsAmqaZAiFAWWCuEOLIS5bpB/QD0jx98TZH6unp8ePHuLm5cfPmTSpWrEijRo1ISEjg0KFDfPTRR6nLJSUlAc8eXuvYsSNRUVE8ffr0nR+Cio2N/c+1Bj8/P1xcXKhRowadO3dO/bxw4cLcunXrnbaTlV2/f52PN32M/xV/GpRpgF8zPyoWqpjh261ZsiZ7uu/hjzN/8PmOz/Fc5MnIOiMZ4zMGc1P1VLSSOaTeoyaESBFCuAElAE9N01xessxCIYSHEMKjUKH/tLM2CM+vAVy7dg0hBHPnzkWv12NtbU1oaGjqKzw8HIDBgwfz6aefcurUKRYsWJDmg1BmZmbo9frU98+Xz5Ur13++e+PGDUxMTLhz585/vpMrV670+pWzhA3hG6g0rxJBN4JY2HIhu7rtypSd/3OaptHRpSPnBp2jR+UeTAyciNdSLyLiIzItg5K9GcRNykKIe8BeoKnsLO/D0tKSH3/8kZkzZ2JpaUmZMmVYs2YN8KztwcmTJwG4f/8+xYsXB2D58uVprtfOzo7jx48DcPz4ca5cuQKAjY0NKSkpqUVAp9PRq1cvfvvtNypWrMisWbNS13HhwgVcXP5TX7OlFH0KI3ePpO0fbXEs4MjJASfp695X2kNpVhZWLGm9hD/a/8H5uPN4/OzBjks7pGRRshdpBUDTtEKapln//edcQCPgnKw86aVKlSq4urry22+/sXLlShYvXkzlypVxdnZm06ZNwLOLuh999BHu7u4ULFgwzXW2a9eO+Ph4nJ2dmTNnDo6Ojqk/a9y4MQcOHABg0qRJeHl5UadOHWbNmsWiRYtSRx179+6lRYsWGfAbG5eHSQ9psaoFUw5OoV/VfgR+HIhDfgfZsQD4yPkjjvc7TimrUjRf1ZwZh2YghFFfFlMMnRBCygtwBU4AYcBpwDet77i7u4t/O3v27H8+y05CQkLE//73v9cu8+TJE1G9enWRnJycSaneTmb9G0Y9jBJVF1QVpuNMxYLgBZmyzXfxMOmhaP9He8FYxKCtg0SKPkV2JMXIAcHiJftUaReBhRBhQBVZ288qqlatSr169UhJSXnlswCRkZFMmTIFM7Ps+9zfxbiLNPm1CXce3eHPzn/SvFxz2ZFeKU+OPKxuv5qvdn3FzMMziU2MZUWbFeQwzSE7mpLFZN89QhbSq1ev1/68XLlylCtXLpPSGJ6zMWept7weeqFnb4+9eBb3lB0pTSaaCTMaz6BI7iJ8tfsr7j65y4aOG7A0t5QdTclCDOIisKJklOc7f1PNlAMfHzCKnf+Lvqz9JUs+WMLuy7tp/XtrnuhU62wl/agCoGRZ4THh1F9eHxPNhD099lC+YHnZkd7Jx1U+Zmnrpfhf9qft6rYk6ZJkR1KyCFUAlCzp6r2rNFjRAIA93fdQoWAFyYneT/fK3VnYaiHbL22nw9oOJKcky46kZAGqAChZTlxiHE1/bcpj3WN2d9+dqQ93ZaQ+Vfswp9kc/jz/J/229FO3iCrvTRWAdGBqaoqbmxsuLi589NFHJCYmvvO6evbsydq1awHo06cPZ8+efeWy+/bt49ChQ2+9DTs7O2JjY9/qO7Nnz2bFihWvXaZTp05cvHjxrfOkp8TkRFr+1pKr967yZ6c/cSmctR5+G+Q5iDE+Y1gWuoxx+8fJjqMYOVUA0sHzVhCnT58mR44czJ8//x8/1+nSnjzkZRYtWoSTk9Mrf/6uBeBt6XQ6lixZQpcuXV673CeffPLOje3SQ4o+hc7rOnPkxhFWtVuFV2kvaVky0hifMXzs9jHj9o9jyYklsuMoRixr3QZqAP2gvby8CAsLY9++fYwePRobGxvOnTtHeHg4I0aMYN++fSQlJTFo0CD69++PEILBgweza9cuSpYsSY4c/3+vd926dZkxYwYeHh789ddfjBo1ipSUFAoWLMjixYuZP38+pqam/Prrr/j5+VGhQgUGDBhAZGQk8OyovXbt2sTFxdG5c2du3rxJzZo1X3nqIE+ePCQkPGvQunbtWrZs2cKyZcvYs2cPVatWxczMDJ1OR82aNZk+fTp169Zl5MiRmJiYMHHiRLy8vOjZsyc6nU7KMwff7PmGP8//iV8zP9pWbJvp288smqaxoOUCbj68Sb/N/ShlVYqG9g1lx1KMUNYqAJLpdDq2b99O06bPWhodP36c06dPU6ZMGRYuXIiVlRXHjh0jKSmJ2rVr07hxY06cOMH58+c5e/Ysd+7cwcnJ6T/39cfExNC3b18CAgIoU6YM8fHx5M+fnwEDBvxjnoAuXbowdOhQ6tSpQ2RkJE2aNCE8PJxx48ZRp04dfH192bp1K4sXL36r3+vgwYO4u7sDzxrTLVu2jPbt2+Pn58dff/3FkSPPmriamJhQtmxZTp48mbp8ZlkZtpKpB6cywH0An3p+mqnblsHc1Jy1H62l1pJadFzbkWN9j2FvYy87lmJkslYBkNQP+nk7aHg2AujduzeHDh3C09MztdXzzp07CQsLSz2/f//+fS5evEhAQACdO3fG1NSUYsWKUb9+/f+sPygoCG9v79R15c+f/6U5du/e/Y9rBg8ePCAhIYGAgADWr18PQIsWLbCxsXmr3y8qKoqKFf//QqqzszPdunWjZcuWHD58+B+jludtpzOzAATfCqbP5j74lPbhh2Y/ZNp2ZcubMy8bO27E42cPPvz9Qw71PkSeHHlkx1KMSNYqAJK8OCXki3Ln/v/p/oQQ+Pn50aRJk38ss23btnTLodfrCQoKwsLi3eabfbEb5ottpl/WdvrUqVNYW1sTHR39j88zu+109KNoPvz9Q4rkLsKaj9Zku3YJDvkdWN1+Nc1WNuPjTR/zR/s/pHU1VYyPugicSZo0acK8efNITn52//aFCxd49OgR3t7erF69mpSUFKKioti7d+9/vlujRg0CAgJS20DHx8cDkDdvXh4+fJi6XOPGjfHz80t9/7woeXt7p04buX37du7evfvSjEWKFCE8PBy9Xs+GDRtSP69YsSKXLl1Kfb9+/Xri4+MJCAhg8ODB3Lt3L/Vnmdl2Wi/0dNvQjdjEWDZ22kih3IY5X0RGa+zQmKkNp7L27FpmHp4pO45iRFQByCR9+vTBycmJqlWr4uLiQv/+/dHpdLRp04Zy5crh5ORE9+7dqVmz5n++W6hQIRYuXEjbtm2pXLkyHTt2BKBVq1Zs2LABNzc3AgMD+fHHHwkODsbV1RUnJ6fUu5HGjBlDQEAAzs7OrF+//pUzq02ZMoWWLVtSq1YtbG1tUz9v1qwZAQEBwLNZyEaMGMGiRYtwdHTk008/5bPPPgPgzp075MqVK3Uayow2OXAyOyN28mOzH3Er6pYp2zRUw2oOo23Ftoz0H8mRG/+ZWE9RXkozpodJPDw8RHBw8D8+Cw8P/8f5aSVjtGnThmnTpr22qdz3339Pvnz56N2791ut+13+Dfdf3U/9FfXp6NyRlW1XqtMewL0n93Cb74amaZzofwJrC+u0v6RkC5qmhQghPP79uRoBKG9kypQpREVFvXYZa2trevTokeFZoh9F03ldZxxsHFjQcoHa+f/N2sKa39v/zo0HN+jzZx/1pLCSJlUAlDdSvnx5vL29X7vMxx9/nOH3/wsh6LGxB/GP41nz0Rry5sybodszNjVK1GBS/UmsC1/HvOB5suMoBk4VAMWozA+ez1+X/mJG4xlULlpZdhyDNKzWMJqVbcYXO77gTPQZ2XEUA6YKgGI0LsRdYNjOYTR2aMygaoNkxzFYJpoJyz5cRt6ceem+sbvqHKq8kioAilHQ6XV029ANCzMLlnywRJ33T0Ph3IVZ0HIBx6OOMzFwouw4ioFSBSAd3L59m06dOuHg4IC7uzvNmzfnwoULUrJMmjTprb+zbNkyPv307don2NnZpf45KiqKli1bvnb5LVu24Ovr+9bZnpsUOImjN48yr8U8iucr/s7ryU7aVmzL/1z/x4SACQTfCk77C0q2owrAexJC0KZNG+rWrUtERAQhISFMnjyZO3fuSMnzLgXgfc2aNYu+ffu+dpkWLVqwefPmd2qVHXIrhPH7x9PZpTMdXTq+a8xsya+ZH0XzFKXbhm48Tn4sO45iYFQBeE979+7F3NycAQMGpH5WuXJlvLy8EELw5Zdf4uLiQqVKlVi9ejXwrI2zj48PrVu3xt7enhEjRrBy5Uo8PT2pVKkSERERwLO5AQYMGICHhweOjo5s2bIF+O8Re8uWLdm3bx8jRoxI7UvUtWtXAH799Vc8PT1xc3Ojf//+pKSkALB06VIcHR3x9PTk4MGDL/3dxo4dy4wZM1Lfu7i4cPXqVeDZw2nPrVu3LrUB3vfff5/azO7UqVO4uLiQmJiIpmnUrVs39Xd4U8kpyfT6sxeFcxdmbvO5b/Vd5dmtoUtbL+Vc7Dm+3fOt7DiKgZHWC0jTtJLACqAIIICFQoj36uT1+V+fE3o7fdtBuxV1Y3bTVzeZO3369Csbn61fv57Q0FBOnjxJbGws1apVS72V8uTJk4SHh5M/f37s7e3p06cPR48e5YcffsDPz4/Zfze2u3r1KkePHiUiIoJ69er9oyXDv02ZMoU5c+aktoAIDw9n9erVHDx4EHNzcwYOHMjKlStp1KgRY8aMISQkBCsrK+rVq0eVKlXe6u/l2LFjAFy5cgUbGxty5swJwGeffUbdunXZsGEDEydOZMGCBVhaWgLg4eFBYGAgHTp0eOPtTDs4jbA7YWzouAGbXG/XxE55ppFDI/q792f2kdl0culEteLVZEdSDITMEYAOGCaEcAJqAIM0TXv17CdG6MCBA6mdPosUKYKPj0/qjrNatWrY2tqSM2dOHBwcaNy4MQCVKlVKPcoG6NChAyYmJpQrVw57e3vOnTv3xtv39/cnJCSEatWq4ebmhr+/P5cvX+bIkSPUrVuXQoUKkSNHjtTWEu8iKirqH6MBExMTli1bRrdu3fDx8aF27dqpP3veKfRNnYs9x/iA8bR3as+HFT5854wKTG04laJ5itJ3c191V5CSStoIQAgRBUT9/eeHmqaFA8WBV8+BmIbXHalnFGdn59QWz2/j+REzPNtpPn9vYmLyjxnE/n23i6ZpmJmZodfrUz/7d6fO54QQ9OjRg8mTJ//j840bN75RxjfZzss6hV68eJE8efL8Z2f/Np1C9UJPnz/7kNs8N37N/NL+gvJaVhZWzG0+lzar2zDz8ExG1BkhO5JiAAziGoCmaXZAFeA/Xaw0TeunaVqwpmnBMTExmR0tTfXr1ycpKYmFCxemfhYWFkZgYCBeXl6pnT5jYmIICAjA09Pzrda/Zs0a9Ho9ERERXL58mfLly2NnZ0doaCh6vZ7r169z9OjR1OXNzc1TO442aNCAtWvXprZsjo+P59q1a1SvXp39+/cTFxdHcnIya9aseem27ezsOH78OPBscpvn3Uhf5Ojo+I8Ry/379xkyZAgBAQHExcX9ozi+TafQecfmcfD6QWY1mUXRPJnTXC6r+7DCh7Sr2I6x+8ZyMU7u3M2KYZBeADRNywOsAz4XQjz498+FEAuFEB5CCI8XTzUYCk3T2LBhA7t378bBwQFnZ2dGjhxJ0aJFadOmDa6urlSuXJn69eszbdq0t+6UWapUKTw9PWnWrBnz58/HwsKC2rVrU6ZMGZycnBgyZAhVq1ZNXb5fv364urrStWtXnJycmDBhAo0bN8bV1ZVGjRoRFRWFra0tY8eOpWbNmtSuXfuVjdjatWtHfHw8zs7OzJkzB0dHx/8skzt3bhwcHFKvTQwdOpRBgwbh6OjI4sWLGTFiRGoB2rt3Ly1atEjzd7754CYj/EfQyL4RPSpnfG+h7MSvmR8WZhb029JP9QpSnp0mkPUCzIEdwBdvsry7u7v4t7Nnz/7ns6yiR48eYs2aNbJjpGn9+vXim2++ee0yt2/fFvXr13/pz/79b9hhTQeR87uc4lLcpXTLqPy/hcELBWMRS08slR1FySRAsHjJPlXaCEB7dnJ7MRAuhJglK4fy/tq0afOPB8NeJjIykpkz056sZPfl3fxx5g9G1hmJQ36HdEqovKh31d7UKlmLr3Z9xd3HL58cSMkepM0HoGlaHSAQOAU8v9I4SgjxyjkS1XwAWdPzf8MkXRKu811J0adweuBpLMzebWpLJW2ht0NxX+jOQI+B+DVXF9mzulfNByDzLqADQLo0dBFCqN4wRurFA5CZh2dyIe4C27tuVzv/DOZW1I2BHgP5KfgnelXpRRXbt3sORMkapF8Efl8WFhbExcWpC1pGSAhBXFwcFhYWXL13lQkBE2hbsS1NyzaVHS1b+K7+dxTIVYBB2wahF/q0v6BkOdJGAOmlRIkS3LhxA0O8RVRJm4WFBSVKlOCjdR+haRqzm2T+sxzZlbWFNdMaTePjTR+zPHQ5H1f5WHYkJZMZfQEwNzenTJkysmMo72H7xe1sOr+JKQ2mUNKqpOw42Ur3yt35+fjPfL37az6s8KFqt5HNGP0pIMW4Jack88XOLyiXvxxDaw6VHSfbMdFMmNt8LnGP4/Dd++7tuhXjpAqAItX84Pmciz3HjMYzyGGaQ3acbMmtqBv9qvZjXvA8wmPCZcdRMpEqAIo08Y/jGbNvDA3KNKCVYyvZcbK18fXGkztHbobvGi47ipKJVAFQpBm/fzz3k+4zq8ksdRuvZIVyF2K092i2XdzGzoidsuMomUQVAEWKc7HnmHtsLn2r9sW1iKvsOAow2HMwDjYOfLHjC3R6XdpfUIyeKgCKFMN3DsfS3JLx9cbLjqL8LadZTqY1msaZmDP8HPKz7DhKJlAFQMl0Oy7tYOvFrYz2Hk3h3IVlx1Fe0KZCG3xK++C7z5d7T+7JjqNkMFUAlEyl0+v4YucXONg4MNhzsOw4yr9omsasJrOIS4xjYsBE2XGUDKYKgJKplp5YytmYs0xvNJ2cZjnT/oKS6araVqWHWw9+PPoj1+5dkx1HyUCqACiZJjE5kTH7xlCrZC01x6+BG193PCaaCb771MNhWZkqAEqm+SHoB6ISopjacKq67dPAlbQqyRDPIfxy8hdO3j4pO46SQVQBUDJFXGIcUw5O4YPyH1CnVB3ZcZQ3MKLOCKwtrBnpP1J2FCWDqAKgZIpJgZNIeJrApPqTZEdR3pBNLhtGeY1i+6Xt7L2yV3YcJQOoAqBkuGv3rjHn2Bx6Vu6Jc2Fn2XGUt/Cp56eUzFeSr3d/rebcyIJUAVAynO8+X0w0E8bWHSs7ivKWLMwsGF9vPMduHWNd+DrZcZR0pgqAkqHC7oTxy8lfGOI5RPX6N1LdXLvhUtiFUf6jSE5Jlh1HSUeqACgZaqT/SKwsrBhRZ4TsKMo7MjUxZUqDKVyMv8ii44tkx1HSkSoASobZd3Uf2y5uY1SdUWqmKSPXvFxzvEt7Mz5gPInJibLjKOlEFQAlQwghGOk/khL5SjC4umr5YOw0TWNCvQncTrjNvGPzZMdR0okqAEqG2HZxG0E3gvD19sXCzEJ2HCUdeJX2oolDE6YcnMLDpIey4yjpQGoB0DRtiaZp0ZqmnZaZQ0lfQghG7x2NvY09Pd16yo6jpKPv6n1HbGIsPxz5QXYUJR3IHgEsA5pKzqCksw3nNnDi9gnG+ozF3NRcdhwlHVUrXo3W5Vsz49AM7j6+++Zf1OshKQkePoT4ePS3biPuRD97r1OTz8hiJnPjQogATdPsZGZQ0leKPgXfvb5UKFiBLpW6yI6jZIDv6n1H5fmVmXl4JhPqTwCdjvuhV7hzKIKEkxHozkdgdusaFg9iyJ0YQ96nsVinxGHC/z9I9u8jTx2mPDLJy/2cRUjIU4QkqyKkFLbFtEJZrDwcsfVxJFf5UmBqmrm/bBYntQC8CU3T+gH9AEqVKiU5jZKW1WdWcybmDKvbr8bURP3PmtWIZB35g5JolOjOrL1TadF1K263z2HFE6z+XiaRXFw3teN+zsJE53YhqWhBdFYFELksIUcOyJkDzdwc9HrE4yfw+DE8eYJpwn1y3r9D3oe3sY47ie2lbeQ59AiWPFtvEjm4buXC/bIemNf0oFgrdwrWqwTmapT5rjTZj3f/PQLYIoRwSWtZDw8PERwcnOGZlHej0+uoOLciluaWnOh/AhNN9hlG5X0JXQrX1x8j6vd9WATtx/72QfKKh5wvAE6DoFNISbpd+Qjh7ELuymUpWN2BUtVtyZP3/bu9JjwUXD1yh+gDF0g8eRHt/DnyR4ZS4VEwNjybreyRlpvLxeqQXKc+JbrVo3DTqmqU8BKapoUIITz+/bnBjwAU47Hi5AouxV9iU6dNaudvxJKi73P+xx08WbeFshe2UUofRynggmlFDtt3BW8fbNvUoOvTcazJ8TvTfxtGsbzF0j1HnrwaLg2LQsOigHfq5wkPBcHbLnNn6zFMDh7A/uoeKq3+GlbDPRMbLpdvjmWnDyg3uCmmNvnSPVdWokYASrpI0iXhOMeRIrmLcKTPEdXv38g8jXvImUmb4PffcLm1E3N0xJGfk8Wbo2vSEvve9XCoWZgX/1mv3L2C4xxH+lXtx9wWc6Vl1+vh7J7bXFu+D3P/v6gatYWCxPEUc84Xr4/WuRMVR7XN1sXgVSMAhBDSXsBvQBSQDNwAer9ueXd3d6EYprlH5wrGIv66+JfsKMobSnnyVJyauEkcte8oHpFLCBCRJqXE9kpfigNTD4hHD3RprmPA5gHCfLy5uHL3SsYHfkP3YpPF7jEBYnP5YeKKVkYIEIlYiOPlO4qIHzYLfdJT2REzHRAsXrJPlT4CeBtqBGCYHic/pqxfWext7AnoGaCO/g1cbPBVLnz1M2UDllA45TaxFOREuQ7k7tsFj8E1yWHx5qfvbjy4Qdkfy9KlUheWtF6SganfzeNEweHvg3iy6Fc8r66mIHHEmRbiar1elJvRn3yVy8iOmCkMcgTwti81AjBMsw7NEoxF7LuyT3YU5RX0yToRNmGTOG7bTKSgCR0m4kD+VmLvsM0i4e77HRF/vv1zYTrOVFyKu5ROaTNG7K0kse2TP8Ve6w+FDhORgibCSjYXET9sFkKX9mjHmPGKEYD0nfrbvFQBMDyPnj4ShacXFg2WN5AdRXmJpLuPRFD3ueJaDgchQNzUiontnr7ign9kum3j1oNbwmKChei1sVe6rTOjndx6XWyq4ituYfvs7yWnnQjt/aPQ3U+QHS1DvKoAqFs1lPcyP3g+0Y+i1WQvBubBxTscauhLQoFSVF8xiHsmBfAf8Ac296/R9Mg4ytVPv7kZbPPa0t+9P8tPLufy3cvptt6M5Nq8BB8cH4fFnWts6bGGWxSn8uIhPLApxfGWvjyJjJYdMXO8rCoY6kuNAAzLo6ePRJHpRUT95fVlR1H+Fn/mljjg+blIxEKkoIkDBVuLw9MDRYpOn6HbfT4K6L2pd4ZuJ6MkJwuxe/xBscf6Q/H8onFwjYEi4dx12dHSBeY+n9kAACAASURBVGoEoKS3hSELufPoDmN8xsiOku3Fn4niYPWhWDjbU/2oH0GlO3FufTi1YzZSY3gdTEwz9sK8bV5b+lXtZ1SjgBeZmUGD0bWoG7+BoKXh7CvelUpBP2NWwYGQ2oNJvHRLdsQMoQqA8k4eJz9m6sGp1LWri3dp77S/oGSIu+G3OVh9KLlc/t7x23Xm8rbz1Lu6FKc25TM1y9d1vsZUM2VS4KRM3W560jSo0bMCzW4s4sz6C+wp3p3Kh+ahlXMgxHsoiZdvy46YrlQBUN7Jz8d/5nbCbXX0L8nj6IccbDQWc6ey/9zxX1mCYzMHKZmK5S1GP/dno4Ard69IyZCeqrSxo9mNnzn5xwX223amcqAfmoM9IU1GkRx7X3a8dKEKgPLWnuieMPXgVLxLe1PXrq7sONlKypNkDnefR4JtWWrvHscJ2+ZEbA6XuuN/0de1v8ZEMzHqUcC/uX9kT9NbSzi+MpwDhdrgvnMyD4uWJazfHMTTZNnx3osqAMpbW3R8Ebce3lJH/5lI6AXHR2/gupULNX8ZyA3L8oTMDcLr1h+Ub1lOdrxUxfMVp1/Vfiw7uYyr967KjpOuPLuUo+GdleyfcYyLOVxw/XkwN6yduTh1PQjjeaD2RaoAKG8lSZfElANTqFOqDvXs6smOky1ErAvlVAEfqk5oS7IwJfCrP3G7vx/3gdVlR3upr+tkvVHAc5oGPsM8qHpvD9s+2UziUzPKjWjHucJe3Np0THa8t6YKgPJWlpxYws2HNxnjM0a1fMhgD67EcaDyQOzau1Psfjj+HRZQ+l4YXlNboZkY7t99iXwl6Fu1L0tDl3Lt3jXZcTKEeQ6N5j+1pFhMGJtaLMAm9hLFPvTkeNU+PI6MkR3vjakCoLyxJF0Skw9MplbJWjQo00B2nCxLn5zCof/9RIpDOWqELWS/y6dw4SINVvcjh6VxdHAfUWdElh0FvCivjRmtt/Qj+cwFtpQfRqUTy3laxpGTffwQyYY/1aUqAMobWxa6jOsPrquj/wx0Zl4Al6zcqbVyEJfzuXFhdSj1T/1AwbLWsqO9lRL5StCnSp8sPQp4UQmnfLQ8N4PjS8M4m8uDyouHcCV/Va6t2C872mupAqC8kacpT5l0YBI1StSgkX0j2XGynPuX4zjk2APngT5YPr1L4JA1VI33x6lDmtNkGKwRdUYAMPnAZMlJMk/1nhWpdncnW3utw+zRA0r3qMuJip15dClKdrSXUgVAeSPLQ5cTeT9SHf2nM6EXHPlsFbqyFah2cRX+1UdhfSscrx/aG/R5/jdR0qokfar2YcmJJUTej5QdJ9OYmWu0WNwWi8tn+dPNl4rnNpDiWIGwT+Y9m73GgKgCoKQpOSWZSQcm4VnckyYOTWTHyTKigq5x3LY51X/sSlQuey6tPk6DoInkKWwpO1q6SR0FBGafUcBzhe0s+eDEOM7+foozltVwnT+Q84Vqc2dXmOxoqVQBUNK04uQKrt67qo7+04k+OYXA9rPJW9OZ8tGB7Gn9AxXiD1GxQyXZ0dJdKatS9KrSi8UnFnP9/nXZcaSo2rEc7nG72NzhF/LHX6JA46qENPyalAePZEdTBUB5veSUZCYGTsSjmAfNyjaTHcfoRWwI41z+mnitG8qZgj7EB5yh/sYhmOU0lR0tw4ysMxKAKQemSE4iT46cGq1W/49HwefYVbwn7v7TuFPIhYg526XmUgVAea1fw37lyr0r6uj/PSXdf0KA1zeUautOoUdXCfxkFZ53tlDKq7TsaBmutHVperr1ZNGJRdx8cFN2HKns3AvQ9Poi/H33k5BigcPg5oSW7yjtIrEqAMor6fQ6JgZOxN3WnRblWsiOY7TC/PYTVcgV7wOTCLLvihYejtdPnY3+Iu/bGOU1Cr3QM/XgVNlRpNM0aDDOm0I3QtlU7TsqXNiErrwTp4cuzvSWEm9UADRNK6xpWhtN0wZpmtZL0zRPTdNU8cjiVoatJOJuBL4+vuro/x3cv3aXAxX74jqkLiZCR/CknXhFLKNg+QKyo2U6O2s7elTuwcKQhdx6mDV7678tm6I5aX30W86sCuN8zsq4zO7DmWINuRsckWkZXrsT1zStnqZpO4CtQDPAFnACvgVOaZo2TtO0fBkfU8lsOr2OCYETqFK0Cq0cW8mOY1yE4MiXa3li70TNc0vY6/ElBW6dxmNk9n5+YpTXKHR6HdMOTpMdxaC4d3bENXYPm5rNp8TtYHJWq0Rot5mZ8yTxy6YJe/4CpgOlXvEzM+BDoN3r1pHG+psC54FLwIi0lldTQmaeX07+IhiLWH92vewoRiUq+IYIsm0tBIizuaqIM7+EyI5kUHpu7CksJliIqIdRsqMYpDM7b4h91h8IAeKClYe4vfNkuqyXV0wJKW1+X8AUiADsgRzAScDpdd9RBSBz6FJ0orxfeeE6z1Wk6FNkxzEKKckpYn+nn8R98opH5BJ7m08TTxOTZccyOBfjLgrTcabii7++kB3FYCU/1Ys//7da3KaweIqZCG76jUh59Pi91vmqAvCm1wB+0TTN6oX3dpqm+b/n4MMTuCSEuCyEeAr8DrR+z3W+VOD0INa2XMrjROPs2Z3Z1pxdw/m484z2Ho2JutSTpojNZzldwBvv3wdy0caT6N2nqLv1S8xzGUfjtsxUNn9Zurp2ZV7wPO4k3JEdxyCZmWu0+qUDicfOsrdoF9z/msj1glWI3BCS7tt60/+7DwBHNE1rrmlaX2AnMPs9t10cePHJkBt/f/YPmqb10zQtWNO04JiYd2uzmnPlYtpv7UVowYYErcy8CyzGSC/0fBfwHc6FnGlbsa3sOAbt6cMk9tUbR8kP3Cj58CwHei+lauwu7BrIn5nLkH3j9Q1JKUnMODRDdhSDVsajAI1uLeevz//iaZLgblIGPCH+smHBy15AHSAZiAKKvun3XrO+9sCiF953A+a87jvvfAooJUWc/Wy+eKDlE4lYiDXVpor4aDU8f5k/Tv8hGIv4/dTvsqMYtLD5B8WlHBWFAHGgdGcRffqO7EhGpeu6rsJyoqWIToiWHcUoJDx4v1OxvOcpoG7AEqA7sAzYpmla5fesPTeBki+8L/H3Z+nPxISKs/tjduEsEeWa0v7Y19wo5on/9OPGOpNbhnh+9F+hYAXaO7WXHccgPbjxgP2VBuE8oA4WKQkcG7uV2ldXUci5sOxoRuVb7295nPyYmYdnyo5iFHLnzaBTsS+rCv9+ARuBwi+89wROvMl3X7NOM+AyUIb/vwjs/LrvpNdF4Ijp60S0WVGhw0RsKDtc3Lz4KF3Wa+zWn10vGIv49eSvsqMYpKBRm8Qtk+IiBU3scxsiHt56IDuSUeu8trPIPTG3iHkUIztKlsf7jACEEB8KIaJfeH8UeK8JSYUQOuBTYAcQDvwhhDjzPut8U/bD22ITFc7p6n348NIMkhwrsfmz3YbWqTVTCSEYHzCecvnL0dGlo+w4BuX2iSgOl/iI6pNak2BuQ/jiw/ic+IE8tnllRzNq33p/S2JyIrMOz5IdJdtK60GwbzVNy/+ynwkhnmqaVl/TtJbvunEhxDYhhKMQwkEIMfFd1/MuzApaUzloATdX7sMslxmtfmzEDtueXDgcl5kxDMaWC1sIvR3KN17fYGai7l4B0Ov07O+yAIuqFalyczP7G36HXVwIzr0MczJ2Y+NUyImPnD/C76gf8Y/jZcfJltIaAZwCNmua5q9p2nRN077SNM3379tCTwGtgCMZHzPjFO/iQ4nYk5xs+Q0No1diXasi6z/6jadJ2efigBCCcfvHYW9jT1fXrrLjGISLG89wKr83Pr8N4IpNVe7sDMNn17eY584hO1qWMtp7NAlPE/j+8Peyo2RLaRWA9kKI2jw7TXOGZw9vPQB+BTyFEEOFEO92b6YB0XJZUHnzBB7sCeFBfjvaru1CUMGWHN+YPWYx2n5pOyFRIeroH3h89wl76vhSuk0VSiaEc7DvUtxi/SndyFF2tCzJpbAL7Z3a8+PRH7n7+K7sONlOWgXAXdO0YkBX4E9gAbACOAbkyuBsma5APVfKRh/mTJ/vcX+0H8c2Tqz1+ZH7cZnQk0MSIQTj94+ntFVpurl2kx1HqpCZ+4gqUpn6B78j2L4j4uw5ai/sma26dsow2ns0D5IeMDvofR8tUt5WWgVgPuAPVACCX3iF/P3frMfUFOefP0ecOkNkKS/aB3xGZNFq7Bh7OEveMrrr8i6O3DzCKK9RmJuay44jRfSZGPaX7Y378HqYi2RCp+6gVsQvFKhQSHa0bMG1iCttKrThhyM/cO/JPdlxspXXFgAhxI9CiIrAEiGE/QuvMkII+0zKKEUe59I4Xd1GxOQ/KKzF0GRcLbYX6835g7Gyo6Wb5+f+S+YrSY/KPWTHyXS6pBT2dpiHuUt5akWsILDmVxS6cxq3rxrLjpbt+Pr4cj/pPj8e+VF2lGzlTW8D/SSjgxgkTcNhxEcUjA7nZKPhNLq9gkJ1HFnXeAEJD4z/ntE9V/Zw6PohRtQZQU6znLLjZKrQBUe4YONJvTUDuWbjxs2tJ/E6NBWL/FlnQnZj4lbUjdblW/N90Pfcf3JfdpxsQ3X6egOm1nmpvHM6DwNCiS7qSrtdA7hUsAa7p4YY9Wmh8QHjKZ63OL2r9JYdJdNEn41lX7m+uA2oQYGkKI4O/Y3Ksf7YNXeSHS3b8/Xx5d6Te/gd9ZMdJdtQBeAt5PdypsKtvVzw/ZUSIpL6I6qxrXgfTu++LTvaW9t/dT8B1wL4uvbX2eLo/2miDv/2P2Hu7EjtS8sI9BxGnpvn8ZzVSV3kNRBVbavS0rElsw7P4mHSQ9lxsgVVAN6WpuE4rivWt88TVv9zGkctp3SjcqyrNpnbV5/ITvfGxgeMp2ieovSp2kd2lAwl9IJDo7cTae1Kg3WDiLSpzI3NoXgdmUHuoupJXkMzxmcMd5/cZc7RObKjZAuqALwjswJWuPnP4nHwWa45NKBd8CiS7CuwsfNqnjw27PNCByIPsOfKHr6q9RW5zLPc3bypzq4+xbGCTak1oTnmJBMyeiOusXso09JZdjTlFTyKedC8XHNmHp5JwtME2XGyPFUA3lM+93K4XNrIjeX+6PNa8+HvnThtU4e/fA+RkiI73cuN3z+ewrkL09+jv+woGeLWiTvsLd+f8p3cKHfvGAfaz6b4vTO4j2+tTvcYAV9vX+IexzH36FzZUbI8VQDSSYnu9SkTH0L4sEWUTrlM0+9qE2jTir2zTxrUheKgG0HsuryL4TWHY2mete54iblwl92eo7Cqak+dC0s47D4Y08uXqLPmM8wsVQsHY1G9RHWaODRhxuEZPHr6SHacLE0VgPRkakrFGb0pEH+Jk50nU+XRAXyGVmFXoS4E/XpJdjoAxu4bS0HLgnxSLevc2XvvRgK76k7EvHwZGh6bzKkyH3DH/wx1gmeTz+6lvQwVAzfGZwyxibHMC54nO0qWpgpABjDJm5vKq0ZgGXWZsOYjqRO/CY9uFdharC9BKyOkjQgORh5kR8QOvqr1FXly5JETIh09uPOYnS1mk1zKnkb7v+VSMR+ubDxJjcu/UaK+6t1jzGqWrEkj+0ZMPzSdxORE2XGyLFUAMpB5YRvctk7E5HIEp2p/QoOoX6j2P0d2Fv4fAfPOZHohGLNvDIVzF2ZgtYGZu+F0Fhtxn531p/DE1o7G24Zys4ArF1ccxuPmJsq0dpUdT0knY3zGEP0omvnB82VHybJUAcgEFnZFqXLADxFxhbAGX1AnbiPeA13YY9OW3VOCM+Vi8f6r+/G/4s+I2iPInSN3xm8wA9w8Ec3Oat9gXrYUjfeO5GbhqpxbsB+3mN2U61ZDdjwlndUuVZsGZRow7eA0NQrIIKoAZKJc9rZU2T0d85vXOPGBLx4P99JwZDWO5a7Lpp7ruRebcV1Hx+wbQ9E8RRngMSDDtpERhIDQFWHscehL/qqlaRg8mfN2TbiyNoQqt7dToZ+37IhKBvL18eXOozv8HPKz7ChZkioAEuSwLUCVTePIE3uNsO4zKM01Wi9vx4PCDmyoOY3zh9N3dqS9V/ay/9p+RtUZZTT3/Sc90rF3yAZC8tXDrUdlalxeSahLN6L8w/G88gdl2lWVHVHJBN6lvalrV5epB6fyRGc8D1oaC1UAJDK1yYfr8mHYPrrE5ZkbeFTEgTZBX1OyVgl2Fv4f24f78/D++zWdE0Lgu8+X4nmL09e9bzolzzjhWyLYUd2XmHz21PNrS7EnlzncZhr6azeoeWohxeuXlx1RyWRjfMYQlRClRgEZ4WUzxRvqy93d/XUT32cJcfvCxPHqA8R9EyshQFzVSouNlX1F0KoIkZLy9uvbeWmnYCxi7tG56R82nURffih2dF4qgnN7CwEiBU2cKNJEhPquE/qnybLjKQbAe6m3KDazmHic/Fh2FKMEBIuX7FM1YUhPKaXBw8NDBAdnzXlo/k0kPubi9I0k/7yUijd3Y4LghHk1rnm0x/bTdnh0dMDUNI11CEGtJbW4+eAmFwdfNKimb3cuPuDU5C3k3LIO95jtWPKYaznLcavxx1SY2A2bSiVkR1QMiP9lfxr+0hC/Zn586vmp7DhGR9O0ECGEx38+VwXA8CWev86Fcb+R+6+1lLt7DIAwsypcdWpBrg+bUOWTGhQs+t+5fLdf3E7zVc1Z0HIB/dz7ZXbsf9Dr4cymS9xcupO8B7bjcXcnOXlKtKktEZXbUmhIZ8p2rwWaatWg/JcQgrrL63Ih7gIRQyKy3FPsGc2gCoCmaR8BY4GKPJtc/o326tm1ALwo8exVLk5dj8W2dZSNDcIUPffJx3HrBjzwqI91s5o4dXKloK0Znos8iU2M5fyn58lhmrmtEHTJgnM7I7m17jBa4H7KXt5JGf1lAG7msCOyahuKDGxHmS410UzVpSglbQciD+C11ItpDafxZe0vZccxKoZWACoCep5NMj9cFYB3kxJ7l8uL9vBw7Q6KndpB0aeRADzCknkuZfiy/RkG3+xNe6eBlG5akZKOuTDJgH1twgM9EbuvELPvDE9PnCHPuWDKxR3GVkQ9y6Pl5kLx+ugaNMF+QGMKVC+rjvSVd9JsZTOO3jzKlc+ukC9nPtlxjIZBFYDUjWvaPlQBSDdPLl7n2m+HeLDrED1dF/HELJFzc8BcDymYcFUrQ0yu0jywLkVSkZJQoiRmRQpgXtCKHIWsyFXUihz5LFLXp2mgf6rj8Z0HPI25z9PYB+hi7qK/fhOzqOtYxl8nf0IkZZIvYMnj1O/dyGFPlF1NqFmTYu1qUqypK5r5f09RKcrbCrkVgsfPHoz1GcuYumNkxzEaqgBkI+vOrqP9mvYsa7mItim1uLXrDAlBp9EunMMyNhKbh5EUTI7ClHe/xfSeiQ1xliV5aFWSpNKOmLs5U9DHmWINnTDLr47MlIzT7o927IrYxZXPrlDAsoDsOEYh0wuApmm7gaIv+dE3QohNfy+zjzQKgKZp/YB+AKVKlXK/du1aBqTNOnR6HZXmVUJDI+yTMMxMXnHknZxM0tUoEiLjSYy6z+PoBzyNvo/+cRJCkNqnSDMzJUfBfOQslA+LIlbkLmZFvgrF0PIafzM5xTidiT5DpXmVGF5rONMaTZMdxyi8qgBk2LhcCNEwndazEFgIz0YA6bHOrGzFyRWciz3Hug7rXr3zBzA3J2e5UuQsVwp1DKUYE+fCznR17cqco3MYWmMotnltZUcyWur2iyzkie4JY/eNpVqxarSp0EZ2HEXJMGN9xpKsT2Zi4ETZUYyalAKgaVobTdNuADWBrZqm7ZCRI6uZHzyf6w+uM7nBZDR1l42ShTnkd6CXWy8Whizk6r2rsuMYLSkFQAixQQhRQgiRUwhRRAjRREaOrORh0kMmBk6koX1DGtg3kB1HUTLcaJ/RmGgmjN8/XnYUo6VOAWURsw7PIjYxlkn1J8mOoiiZokS+EgysNpDlJ5dzLvac7DhGSRWALCDmUQwzDs+gbcW2VCteTXYcRck0I+uMJJdZLkbvHS07ilFSBSALmHxgMonJiUyoN0F2FEXJVIVyF2J4reGsPbuWIzeOyI5jdFQBMHKR9yOZe2wuPSr3oGKhirLjKEqmG1ZzGEVyF+HLXV9iTM0tDYEqAEZu7L6xz/5bd6zUHIoiS96ceRnjM4bAyEC2XNgiO45RUQXAiIXdCWNZ6DI+rfYppaxKyY6jKNL0qdoHxwKOfL37a3T6jJtbO6tRBcCIfbnrS6wtrPnW+1vZURRFKnNTcyY3mEx4bDjLQpfJjmM0VAEwUn9d+oudETvx9fHFJpeN7DiKIl2bCm2oWaImvnt9efT0kew4RkEVACOk0+sYvnM4DjYODKw2UHYcRTEImqYxvdF0ohKimB00W3Yco6AKgBFaFrqMMzFnmNJwSqbP9KUohqx2qdp8WOFDph6cSsyjGNlxDJ4qAEYm4WkCo/eOplbJWrSr2E52HEUxOJMbPHsuZsw+NWFMWlQBMDLTD07ndsJtZjaeqRq+KcpLVChYgU88PmFByAJO3TklO45BUwXAiNx8cJPph6bTwbkDNUrUkB1HUQzWuHrjsLaw5vMdn6uHw15DFQAj8s2eb9DpdUxuMFl2FEUxaPlz5Wdc3XHsubKHTec3yY5jsFQBMBJBN4JYfnI5X9T8Ansbe9lxFMXgDfAYgFMhJ4bvHE6SLkl2HIOkCoAR0As9g7cPxjaPLd94fSM7jqIYBTMTM2Y3mU3E3Qh+OPKD7DgGSRUAI7D0xFKCbwUzvdF08ubMKzuOohiNRg6NaOXYiu8CvuN2wm3ZcQyOKgAG7t6Te4z0H0ntkrXpUqmL7DiKYnRmNp5Jki6Jb/zV6PnfVAEwcGP3jSU2MRa/Zn7qtk9FeQflCpTjs+qfsTT02Uha+X+qABiwM9FnmHN0Dv3c+1HFtorsOIpitL71/pYieYrwydZPSNGnyI5jMFQBMFBCCIb8NYR8OfMxsf5E2XEUxahZWVgxq/Esgm8FMz94vuw4BkMVAAP1a9iv7Lmyh0kNJlHAsoDsOIpi9Dq5dKKhfUNG7RmlLgj/TRUAAxSXGMcXO7+gRoka9HPvJzuOomQJmqYxt/lcnuieMGznMNlxDIKUAqBp2nRN085pmhamadoGTdOsZeQwVF/t+op7T+6xsOVCTDRVoxUlvTgWcGRE7RGsOrUK/8v+suNIJ2vvsgtwEUK4AheAkZJyGJz9V/ezJHQJw2sOp1KRSrLjKEqWM9Jr5LO5NLYNzPZPCEspAEKInUKI5xN3BgElZOQwNEm6JPpv6U8Z6zKM9hktO46iZEkWZhbMbT6XC3EXmHZwmuw4UhnC+YVewPZX/VDTtH6apgVrmhYcE5O1J3iYenAq5+PO81OLn7A0t5QdR1GyrCZlm9DJpRMTAidwJvqM7DjSZFgB0DRtt6Zpp1/yav3CMt8AOmDlq9YjhFgohPAQQngUKlQoo+JKdzbmLBMDJ9LJpRNNyzaVHUdRsrwfm/5Ivpz56PVnr2z7bECGFQAhREMhhMtLXpsANE3rCbQEuops3rBbp9fRc2NP8ubIy+wmai5TRckMhXIXwq+ZH0dvHs22cwjLuguoKfAV8IEQIlFGBkMy49AMjt06xk8tfqJIniKy4yhKttHRuSOty7fm273fcjHuouw4mU7WNYA5QF5gl6ZpoZqmZdtH805Hn2bMvjG0d2pPB+cOsuMoSraiaRo/tfiJnKY56f1nb/RCLztSppJ1F1BZIURJIYTb368BMnLIlpySTM+NPbHKacVPzX+SHUdRsqVieYvxfZPvCYwM5Ieg7DVvgCHcBZRtTT04lZCoEOa1mEeh3Fn3AreiGLqebj35oPwHjPAfka0mklcFQJIjN44wdt9YOrl0op1TO9lxFCVb0zSNRa0WYWNhQ9f1XXmieyI7UqZQBUCCB0kP6LyuMyXylWBei3my4yiKwrO7gpa2Xsqp6FPZZvIYVQAymRCCT7Z+QuT9SFa1W4W1hWqDpCiGolm5ZgyqNohZQbOyRa8gVQAy2S9hv7Dq1CrG+IyhVslasuMoivIv0xpNo0LBCnTf2J3oR9Gy42QoVQAy0cW4iwzaNgjv0t6M8holO46iKC9haW7J6variX8cT5d1XbL0U8KqAGSSR08f0e6PduQwzcGvbX7F1MRUdiRFUV7BtYgr81rMw/+KP2P3jZUdJ8OoApAJhBD029KP09GnWdV2FSWtSsqOpChKGnq69aR3ld5MCJzA9ouv7Fdp1FQByARzj81l1alVjK83niZlm8iOoyjKG/Jr5kflIpX534b/ce3eNdlx0p0qABnsQOQBhu4YSivHVuq8v6IYmVzmuVjbYS06vY4PV39IwtME2ZHSlSoAGejy3cu0Wd2GMtZlWNFmhZreUVGMUNn8ZVndfjVhd8LotqFbluoXpPZIGeTek3u0WNUCvdCztctWdb+/ohixpmWb8n2T79l4bmOWekjMTHaArCg5JZn2f7QnIj6CXd12Ua5AOdmRFEV5T4M9B3M25ixTDk7BsYAjH1f5WHak96YKQDrTCz19NvfB/4o/y1ovw8fOR3YkRVHSgaZp+DXzI+JuBH0396WgZUFalW8lO9Z7UaeA0pEQgs//+pwVJ1cwvu54erj1kB1JUZR0ZG5qzvoO66liW4UOazsQeC1QdqT3ogpAOhqzbwx+R/0YVnMY33p/KzuOoigZIG/OvGzrso3SVqVp+VtLQm+Hyo70zlQBSCdTD0zlu4Dv6F2lN9MbTUfTNNmRFEXJIIVyF2Jnt53ky5mPhisaGm0RUAXgPQkh8N3rywj/EXR26cyClgvUzl9RsoFSVqXY22MvluaW1F9en+BbwbIjvTVVAN6DEILhO4fzXcB39HLrxS9tflE9fhQlGymbvywBHwdgZWFFgxUNOHz9sOxIb0UVgHeUnJJMnz/7MCtoFoM9B/PzBz+rnb+iZEN21nYE9AygcO7CNFjRgPXh62VHemOqALyD+dh5HwAACoxJREFUu4/v0nRlU5aELmG092h+aPqDespXUbKxklYlOdjrIJWLVqb9H+2ZcWgGQgjZsdKk9lpv6Uz0GWourkngtUCWf7ic8fXGq3P+iqJQOHdh9nTfQ3un9ny560v6bu5r8HMLSykAmqZ9p2lamKZpoZqm7dQ0rZiMHG9reehyqv1cjXtP7rG7+266V+4uO5KiKAYkl3kufm//O6PqjGLxicXUWFSDi3EXZcd6JVkjgOlCCFchhBuwBfCVlOONxCbG0mVdF3pu6kn1EtUJHRCKd2lv2bEURTFAJpoJExtMZEvnLVx/cB33he4sDFlokKeEpBQAIcSDF97mBgzvb4Znd/n8fvp3nOY6sfbsWsbXHc/ubrspmqeo7GiKohi4Fo4tCO0fikcxD/pv6U/DXxoSER8hO9Y/SLsGoGnaRE3TrgNdec0IQNO0fpqmBWuaFhwTE5Np+Y7cOILXUi86r+tMaevShPQLYbTPaHWnj6Iob6ykVUn8u/uzoOUCjt08RsW5FRm2Yxh3H9+VHQ0ALaOGJZqm7QZedqj8jRBi0wvLjQQshBBj0lqnh4eHCA7OuIcthBAEXAtg+qHpbL24lSK5i/Bdve/4uMrHmJmovnmKory7Ww9vMXrPaJaGLsXKwoqBHgMZXH1wppxR0DQtRAjh8Z/PZZ+X0jStFLBNCOGS1rIZVQCiHkbx2+nf+CXsF0Jvh1LIshBDqg/hs+qfkTdn3nTfnqIo/9fe3cVIdZdxHP8+LuwuWF0GlrRgCbANwhJfoMGmwURebGltzIIt6JqQbhUTqdobYyINFzVEU/WmCdGkmAbrS4TWbQhr+kKgu6QXlrW9gEJZYRfQCkJZQSBGWSg8Xpz/NIdlZndm58zMsuf3SSZz5n/efvucs/Ofc+bMTHodOHOATW9sYkfPDmpranm4+WHWfmYt9zfdz/ia8WVZ56jqAMxsjrv3huEngCXuvnq4+ZLoAAY+GKD3fC89/T10n+qm80Qn+8/sx3EWTV/EuoXraPtsGxPGTyhpPSIiQzl67iibuzez7dA2zv/vPJn6DEtnLWXZrGUsuGMBzVObaZzYmMi6RlsH8BIwF7gO/B1Y7+6nhptvpB3AU11PsXX/Vi5cvnDDb3rW1dSxeMZils9ezpr5a5jbOLfoZYuIlOLKtSu81vcaO/+6k66/dXHiwokPx00cP5GGugYa6hvY8uUtI776MF8HUJUT2+7+SCXX15RpYkXTCibVTyIzIcNdmbuY1ziP5qnN1I+rr2QUEZEb1NbU0jK3hZa5LQC8d/E9Dvcfpqe/h5OXTnJp4BIXBy7SUNeQ+Lqr/h5AMcr9JrCIyFiU7whAXwUhIpJS6gBERFJKHYCISEqpAxARSSl1ACIiKaUOQEQkpdQBiIiklDoAEZGUuqU+CGZm/URfHTESjcC/EoyTFOUqjnIVR7mKM1pzQWnZZrr71MGNt1QHUAozezvXJ+GqTbmKo1zFUa7ijNZcUJ5sOgUkIpJS6gBERFIqTR3Ar6odIA/lKo5yFUe5ijNac0EZsqXmPQAREblRmo4AREQkRh2AiEhKjakOwMzWmNm7ZnbdzPJeLmVmD5rZETPrM7MNsfbZZtYd2l8ws9qEck02s91m1hvuMzmmWWZm+2O3y2a2Kox73sxOxMYtqFSuMN212Lo7Yu3VrNcCM3szbO93zOxrsXGJ1ivf/hIbXxf+/r5Qj1mxcU+G9iNm9kApOUaQ6/tmdjjU53Uzmxkbl3ObVijXY2bWH1v/t2Lj2sJ27zWztgrneiaW6aiZXYiNK2e9tprZWTM7lGe8mdnmkPsdM7s7Nq60ern7mLkBzUS/NbwXWJRnmhrgGNAE1AIHgPlh3ItAaxh+Fng8oVw/BzaE4Q3Az4aZfjJwHpgYHj8PrC5DvQrKBfwnT3vV6gV8EpgThqcDp4FJSddrqP0lNs13gGfDcCvwQhieH6avA2aH5dRUMNey2D70eDbXUNu0QrkeA36RY97JwPFwnwnDmUrlGjT9E8DWctcrLPsLwN3AoTzjHwJeBQy4F+hOql5j6gjA3Xvc/cgwk90D9Ln7cXe/AmwHVpqZAcuB9jDdb4BVCUVbGZZX6HJXA6+6+38TWn8+xeb6ULXr5e5H3b03DP8TOAvc9EnHBOTcX4bI2w58MdRnJbDd3Qfc/QTQF5ZXkVzu3hXbh/YBdya07pJyDeEBYLe7n3f3fwO7gQerlOvrwLaE1j0kd3+D6AVfPiuB33pkHzDJzKaRQL3GVAdQoE8A/4g9PhnapgAX3P2DQe1JuN3dT4fhM8Dtw0zfys0730/C4d8zZlZX4Vz1Zva2me3LnpZiFNXLzO4helV3LNacVL3y7S85pwn1uEhUn0LmLWeuuHVEryKzcm3TSuZ6JGyfdjObUeS85cxFOFU2G+iMNZerXoXIl73keo0rOVqFmdke4I4coza6+85K58kaKlf8gbu7meW99jb07J8GdsWanyR6Iqwluhb4h8CmCuaa6e6nzKwJ6DSzg0RPciOWcL1+B7S5+/XQPOJ6jUVmthZYBCyJNd+0Td39WO4lJO5PwDZ3HzCzbxMdPS2v0LoL0Qq0u/u1WFs161U2t1wH4O73lbiIU8CM2OM7Q9s5okOrceFVXLa95Fxm9r6ZTXP30+EJ6+wQi/oqsMPdr8aWnX01PGBmvwZ+UMlc7n4q3B83s73AQuAlqlwvM/s48DJR578vtuwR1yuHfPtLrmlOmtk4oIFofypk3nLmwszuI+pUl7j7QLY9zzZN4glt2Fzufi728Dmi93yy8y4dNO/eBDIVlCumFfhuvKGM9SpEvuwl1yuNp4DeAuZYdAVLLdHG7vDoXZUuovPvAG1AUkcUHWF5hSz3pnOP4Ukwe959FZDzaoFy5DKzTPYUipk1Ap8HDle7XmHb7SA6N9o+aFyS9cq5vwyRdzXQGerTAbRadJXQbGAO8JcSshSVy8wWAluAFnc/G2vPuU0rmGta7GEL0BOGdwErQr4MsIIbj4TLmitkm0f0huqbsbZy1qsQHcCj4Wqge4GL4UVO6fUq1zvb1bgBXyE6DzYAvA/sCu3TgVdi0z0EHCXqwTfG2puI/kH7gD8CdQnlmgK8DvQCe4DJoX0R8FxsullEvfpHBs3fCRwkeiL7PXBbpXIBi8O6D4T7daOhXsBa4CqwP3ZbUI565dpfiE4ptYTh+vD394V6NMXm3RjmOwJ8KeH9fbhce8L/QbY+HcNt0wrlehp4N6y/C5gXm/eboY59wDcqmSs8/hHw00Hzlbte24iuYrtK9Py1DlgPrA/jDfhlyH2Q2BWOpdZLXwUhIpJSaTwFJCIiqAMQEUktdQAiIimlDkBEJKXUAYiIpJQ6ABGRlFIHICKSUuoAREpgZp8LX2pWb2Yftej3CT5V7VwihdAHwURKZGY/Jvo08ATgpLs/XeVIIgVRByBSovDdMm8Bl4HFfuO3SIqMWjoFJFK6KcBtwMeIjgREbgk6AhApkUW/Ebud6EdEprn796ocSaQgt9zvAYiMJmb2KHDV3f9gZjXAn81subt3DjevSLXpCEBEJKX0HoCISEqpAxARSSl1ACIiKaUOQEQkpdQBiIiklDoAEZGUUgcgIpJS/wd4aNChXP8F1QAAAABJRU5ErkJggg==\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        ""
      ],
      "metadata": {
        "id": "gTIEZ311xCVK"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}
